In this experiment, we played amplitude modulated white noise to a group of bats in order to investigate their ability to temporally modulate their calling behavior at the group level.
We hypothesized that bats would adapt the timing of their calls to avoid temporal overlap with the amplitude modulated noise.
After extracting call timins from acoustic recordings (.wav files) and identifying the instantaneous period (or phase) in the real or simulated (as in the silent condition) amplitude modulation cycle, we analyzed the call onset data, primarily with circular statistics.
packages <- c("circular", "CircStats", "CircSpaceTime", "sjPlot", "yarrr", "janitor", "broom", "knitr", "kableExtra", "car", "MASS", "emmeans", "rcompanion", "pscl", "randomForest", "chisq.posthoc.test", "viridis", "scales","ggh4x","ggside", "ggdist", "ggridges","gganimate","plotly","htmlwidgets","tidyverse")
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
str(params)
## List of 7
## $ subsample : logi FALSE
## $ data_path : chr "./../1-data"
## $ results_path: chr "./../3-results"
## $ exp : chr "exp_2"
## $ plot_sm : logi TRUE
## $ use_all : logi TRUE
## $ seed : num 42
# local paths
data_path <- params$data_path
results_path <- params$results_path
tables_path <- file.path(results_path,'tables')
if (!dir.exists(tables_path)) dir.create(tables_path)
figures_path <- file.path(results_path,'figures')
if (!dir.exists(figures_path)) dir.create(figures_path)
models_path <- file.path(results_path,'models')
if (!dir.exists(models_path)) dir.create(models_path)
sm_path <- file.path(results_path,'supplementary_figures')
if (!dir.exists(sm_path)) dir.create(sm_path)
# experiment number
exp <- params$exp
# seed for subsampling
seed <- params$seed
# whether to use all data for circular stats
use_all <- params$use_all
subsample <- params$subsample
# whether to plot supplementary plots/tables
supplementary <- params$plot_sm
# load data
data <- readr::read_rds(file.path(data_path,'input',paste0(exp, "_data.RDS"))) %>%
mutate(start_phase_circ = circular(start_phase, units="radians", zero = 0, modulo = "2pi", rotation = "counter"),.after ="start_phase")
# data for averaging histograms (modulation variable is numeric-ish)
data_2proc <- data %>%
mutate(modulation = factor(modulation,
levels = sort(unique(data$modulation))))
# data for remaining analyses (modulation variable is character-ish, better for plotting)
data <- data %>%
mutate(modulation = factor(modulation,
levels = sort(unique(data$modulation)),
labels = paste0(as.character(sort(unique(data$modulation))),"Hz")))
str(data)
## 'data.frame': 1088994 obs. of 23 variables:
## $ group : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ session : int 1 1 1 1 1 1 1 1 1 1 ...
## $ part : int 1 1 1 1 1 1 1 1 1 1 ...
## $ condition : Factor w/ 3 levels "silence","steady-state masker",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ modulation : Factor w/ 8 levels "4Hz","8Hz","16Hz",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ minute : num 1 1 1 1 1 1 1 1 1 1 ...
## $ start_seconds : num 0.0151 0.0777 0.1069 0.1599 0.1902 ...
## $ stop_seconds : num 0.0169 0.0794 0.1086 0.1627 0.1918 ...
## $ start_sp : num 3776 19437 26717 39986 47554 ...
## $ end_sp : num 4223 19848 27141 40671 47960 ...
## $ start_cycle : int 1 2 2 3 4 4 5 6 7 9 ...
## $ end_cycle : int 1 2 2 3 4 4 5 6 7 9 ...
## $ start_phase : num 1.518 1.533 4.46 3.513 0.273 ...
## $ start_phase_circ: 'circular' num 1.518 1.533 4.46 3.513 0.273 ...
## ..- attr(*, "circularp")=List of 6
## .. ..$ type : chr "angles"
## .. ..$ units : chr "radians"
## .. ..$ template: chr "none"
## .. ..$ modulo : chr "2pi"
## .. ..$ zero : num 0
## .. ..$ rotation: chr "counter"
## $ end_phase : num 1.698 1.698 4.631 3.788 0.436 ...
## $ start_per : num 0.0151 0.01524 0.04437 0.03494 0.00271 ...
## $ end_per : num 0.01689 0.01689 0.04606 0.03768 0.00434 ...
## $ duration : num 0.00179 0.00164 0.0017 0.00274 0.00162 ...
## $ onset_interval : num NA 0.0626 0.0291 0.0531 0.0303 ...
## $ class : chr "call" "call" "call" "call" ...
## $ file_no : int 1 1 1 1 1 1 1 1 1 1 ...
## $ familiarization : num 0 0 0 0 0 0 0 0 0 0 ...
## $ order : num 1 1 1 1 1 1 1 1 1 1 ...
# save interim data manips to 1-data/output
data_path <- file.path(data_path,'output')
if (!dir.exists(data_path)) dir.create(data_path)
A few rough summary statistics. Note that each row is an observation (call event), with parameters:
#summary(data)
Total number of calls:
nrow(data)
## [1] 1088994
#1,425,067 # exp 1
#1,088,994 # exp 2
# = 2,514,061 # total
Call durations:
# durations
# data1 <- data
# data2 <- data
# data <- rbind(data1,data2)
# nrow(data)
summary(data$duration)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000996 0.001716 0.003128 0.003881 0.004364 0.704744
# # Min. 1st Qu. Median Mean 3rd Qu. Max.
# #0.000996 0.001704 0.003372 0.004246 0.005036 0.704744
IQR(data$duration)
## [1] 0.002648
# # 0.003332
Calculate histogram averages (essentially, a PSTH)
call_summary1 <- get_avg_hist(data_2proc, NULL, list("group","session"), "period",
by_group = TRUE, by_session = TRUE)
call_summary60 <- get_avg_hist(data_2proc, 60, list("group","session"), "period")
# optionally: check that the two have the correct # of bins
# call_summary1$generic %>% group_by(modulation) %>% summarise(sort(unique(mids),decreasing = T)[2])
# call_summary60$generic %>% group_by(modulation) %>% summarise(n=length(unique(mids)))
First, we asked: Did the presence of playback noise affect the distribution of emitted calls within the modulation cycle?
Here we plot an averaged histogram of detected call onsets (binned calls (1 ms bins) summed over all cycles, averaged over groups and recording days):
Points are mean # calls, shaded areas are +- sem…
# average call counts all at once, 1 ms bins
p <- p.hist_avg_dot_kwargs(dat = call_summary1$generic, facets = c("condition","modulation"))
print(p)
ggsave(file.path(figures_path,paste0(exp, '_fig1a.png') ), p, width = w, height = 5)
# average call counts with fixed y-axis, to better illustrate the relative calling rate between conditions
p2 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary1$generic, facets = c("","modulation"))
print(p2)
ggsave(file.path(figures_path,paste0(exp, '_fig1a-alt1.png') ), p2, width = w, height = 5)
# average call counts with fixed y-axis, with 60 ms bins for each modulation rate, to show trend on the same x-scale across conditions
p3 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary60$generic, facets = c("","modulation"))
print(p3)
ggsave(file.path(figures_path,paste0(exp, '_fig1a-alt2.png') ), p3, width = w, height = 5)
# show average call counts (with and without fixed y axis) for each group, and for each recording day ("session")
for (m in sort(unique(call_summary1$group$modulation))) {
p1 <- p.hist_avg_dot_kwargs(dat = call_summary1$group %>% filter(modulation == m),
facets = c("condition","group")) +
ggtitle(paste0(m, ": by group"))
p2 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary1$group %>% filter(modulation == m),
facets = c("","group")) +
ggtitle(paste0(m, ": by group"))
p3 <- p.hist_avg_dot_kwargs(dat = call_summary1$session %>% filter(modulation == m),
facets = c("condition","session")) +
ggtitle(paste0(m, ": by day"))
print(p1); print(p2); print(p3)
ggsave(file.path(sm_path,paste0(exp, '_figs3a_',m,'.png') ),
p1, width = 7, height = 5)
ggsave(file.path(sm_path,paste0(exp, '_figs3a_',m,'-alt1.png') ),
p2, width = 7, height = 5)
ggsave(file.path(sm_path,paste0(exp, '_figs3b_',m,'.png') ),
p3, width = 7, height = 5)
}
# show average hist of binned calls for the first 7.5-8 mins of continuous AM playback for each mod.
call_summary1_init <- get_avg_hist(data_init, NULL, list("group"), "period", by_group = TRUE)
#str(data_init)
firsts <- p.hist_avg_col(dat = call_summary1_init$group %>%
dplyr::filter(condition=="full-band masker" | condition=="steady-state masker"),
facets = c("modulation","group"))
print(firsts)
ggsave(file.path(sm_path,paste0(exp, '_figs4.png') ), firsts, width = 7, height = h2)
To analyze the distribution of call onsets within the modulation
cycle, we must use circular statistics (with functions from the package
circular), since the nature of our data is cyclical.
First we will compute the following summary statistics:
circular::mle.vonmises())circular::mle.vonmises())circular::mle.vonmises.bootstrap.ci())Then, we’ll get bootsrapped mean and concentration parameters from the whole dataset.
We predicted that calls emitted in presence of modulated noise would tend to “cluster” unimodally in the hemisphere corresponding to the amplitude trough. Conversely, we predicted that the distribution of call onsets emitted in silence would show no clustering, and therefore resemble a uniform distribution.
This prediction entails two specific hypotheses:
Distribution of calls detected in silence will be equivalent to a
uniform distribution, while that of calls detected in the presence of
masking noise will deviate from the uniform distribution. –> Method:
Rayleigh test of uniformity
(circular::rayleigh.test(x, mu=circular(0)))
The distribution of calls in the playback conditions will each
besignificantly different from the silence condition. –> Method:
Mardia-Watson-Wheeler non-parametric test of homogeneous samples
(circular::watson.wheeler.test(x,y)), where differences
between samples can be in either the mean or the variance.
In which conditions were call onset distributions unimodally clustered?
# prepare data for analysis: split into lists, one for each modulation x condition case
# data is subsampled if requested in header
data_split <- data_samp %>% group_by(modulation, condition) %>% group_split()
# label and convert to circular data
circ_data <- lapply(data_split, function(x)
list(modulation = unique(x$modulation),
condition = unique(x$condition),
phase = x$start_phase_circ))
if (file.exists(f_circsum)) { # as RDS...
circ_summary <- readr::read_rds(f_circsum)
} else {
# calculate values:
circ_datsum <- lapply(circ_data,
function(x) tibble(modulation = x$modulation,
condition = x$condition,
n = length(x$phase),
theta_bar = mean.circular(x$phase, na.rm = T), # angular mean
r_bar = rho.circular(x$phase, na.rm = T), # mean resultant length
vm = var.circular(x$phase),
v = sd.circular(x$phase),
# Test of uniformity:
rt = rayleigh.test(x$phase, mu=circular(0))$statistic,
p = rayleigh.test(x$phase, mu=circular(0))$p.value,
# MLE von Mises:
mle_mu = mle.vonmises(x$phase)$mu,
bs_mu_ci_l = mle.vonmises.bootstrap.ci(x$phase, alpha=.05)$mu.ci[1],
bs_mu_ci_h = mle.vonmises.bootstrap.ci(x$phase, alpha=.05)$mu.ci[2],
mle_kappa = mle.vonmises(x$phase)$kappa)
)
# wrangle & correct p values
circ_summary <- tibble(do.call(rbind,circ_datsum)) %>%
mutate(p_adjust = p.adjust(p, method ="bonferroni"), .after = "p") %>%
mutate(across(where(is.numeric), round, 3))
head(circ_summary)
# save, if first time
if (subsample) {
write.csv(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data_seed",seed,".csv")), row.names = F)
readr::write_rds(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data_seed",seed,".RDS")))
} else {
write.csv(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data-alldata.csv")), row.names = F)
readr::write_rds(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data-alldata.RDS")))
}
}
circ_summary
## # A tibble: 24 × 14
## modulation condi…¹ n theta…² r_bar vm v rt p p_adj…³ mle_mu bs_mu…⁴ bs_mu…⁵ mle_k…⁶
## <fct> <fct> <dbl> <circu> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <circ> <circu> <circu> <dbl>
## 1 4Hz silence 59635 6.203 0.003 0.997 3.39 0.003 0.134 1 6.203 4.207 1.704 0.006
## 2 4Hz steady… 28938 5.246 0.074 0.926 2.28 0.038 0 0 5.246 5.132 5.353 0.148
## 3 4Hz random… 105990 5.245 0.054 0.946 2.42 0.027 0 0 5.245 5.165 5.326 0.108
## 4 8Hz silence 62203 5.480 0.002 0.998 3.46 0.002 0.27 1 5.480 2.943 2.224 0.005
## 5 8Hz steady… 65446 5.795 0.057 0.943 2.40 0.05 0 0 5.795 5.699 5.890 0.114
## 6 8Hz random… 51939 5.459 0.061 0.939 2.37 0.041 0 0 5.459 5.366 5.567 0.122
## 7 16Hz silence 65769 0.858 0.003 0.997 3.38 0.002 0.213 1 0.858 5.194 3.164 0.007
## 8 16Hz steady… 60482 5.732 0.042 0.958 2.52 0.036 0 0 5.732 5.610 5.851 0.084
## 9 16Hz random… 25543 5.376 0.046 0.954 2.49 0.028 0 0 5.376 5.178 5.571 0.091
## 10 25Hz silence 67403 4.819 0.003 0.997 3.37 0 0.448 1 4.819 2.495 1.907 0.007
## # … with 14 more rows, and abbreviated variable names ¹condition, ²theta_bar, ³p_adjust, ⁴bs_mu_ci_l,
## # ⁵bs_mu_ci_h, ⁶mle_kappa
Now, let’s plot our call count data (as a density rather than raw call count) as before, but on the polar plane and indicate the mean direction and resultant length of the resultant vector.
circ_density <- p.circ_density(dat = data_samp, circ_dat = circ_summary_plot, bins = 30) %>% recolor(which="f")
circ_density
if (subsample) {
ggsave(file.path(figures_path,paste0(exp, "_fig1b_",seed,".png")), circ_density, width = w, height = 5)
} else {
ggsave(file.path(figures_path,paste0(exp, '_fig1b.png')), circ_density, width = w, height = 5)
}
How consistent was the clustering, if present?
circ_params <- p.circ_params(mu_bs_df, compress = T) %>% recolor()
circ_params
ggsave(file.path(figures_path,paste0(exp, '_fig1c.png')), circ_params, width = 3, height = 6)
mu_kappa <- p.mu_kappa(dat = mu_bs_df %>%
group_by(modulation, condition) %>%
slice_sample(n = 500)) %>% recolor()
mu_kappa
ggsave(file.path(figures_path,paste0(exp, '_fig1d-500.png') ), mu_kappa, width = 6, height = 6)
We now compute frequentist non-parametric analyses for whether the distributions of calls observed in our different modulation and masking conditions are identical.
The null hypothesis is that the distributions are identical, either in their means or their variances:
Did distributions of call onsets differ between playback conditions at each modulation rate?
circ_data_df <- do.call(rbind, lapply(circ_data, function(x) as.data.frame(x)))
ww_tests <- lapply(sort(unique(circ_data_df$modulation)), function(x)
watson.wheeler.test(phase ~ condition,circ_data_df[circ_data_df$modulation==x,]) %>%
broom::tidy() %>%
mutate(across(where(is.numeric), round, 3)) %>%
mutate(modulation = x,.before=1)) %>%
do.call(rbind,.) %>%
mutate(p.adjust = p.adjust(p.value, method ="bonferroni"), .after = "p.value") %>%
janitor::clean_names()
ww_tests
## # A tibble: 8 × 6
## modulation statistic p_value p_adjust parameter method
## <fct> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 4Hz 282. 0 0 4 Watson-Wheeler test for homogeneity of angles
## 2 8Hz 274. 0 0 4 Watson-Wheeler test for homogeneity of angles
## 3 16Hz 143. 0 0 4 Watson-Wheeler test for homogeneity of angles
## 4 25Hz 45.4 0 0 4 Watson-Wheeler test for homogeneity of angles
## 5 33Hz 32.7 0 0 4 Watson-Wheeler test for homogeneity of angles
## 6 40Hz 20.6 0 0 4 Watson-Wheeler test for homogeneity of angles
## 7 50Hz 8.78 0.067 0.536 4 Watson-Wheeler test for homogeneity of angles
## 8 80Hz 2.29 0.683 1 4 Watson-Wheeler test for homogeneity of angles
Test if means or vector lengths were significantly different between conditions and, separately, across modulation rates:
How did distributions of call onsets differ between playback conditions at each modulation rate? In their angular means or their polar concentrations?
rao_testsm <- do.call(rbind, lapply(sort(unique(circ_data_df$modulation)), function(x) {
df <- circ_data_df[circ_data_df$modulation==x,] %>% group_split(condition)
rao <- circular::rao.test(lapply(df,'[[', 'phase'))
data.frame(modulation = as.factor(x),
test = as.factor(c("polar vectors", "dispersions")),
statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>%
mutate(across(where(is.numeric), round, 3)) %>%
mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value")
}))
# test differences for playback conditions within mod rates
rao_testsm
## modulation test statistic df p_value p_adjust
## 1 4Hz polar vectors 3.110 2 0.211 0.422
## 2 4Hz dispersions 232.370 2 0.000 0.000
## 3 8Hz polar vectors 18.883 2 0.000 0.000
## 4 8Hz dispersions 199.932 2 0.000 0.000
## 5 16Hz polar vectors 6.895 2 0.032 0.064
## 6 16Hz dispersions 77.400 2 0.000 0.000
## 7 25Hz polar vectors 3.853 2 0.146 0.292
## 8 25Hz dispersions 22.588 2 0.000 0.000
## 9 33Hz polar vectors 1.668 2 0.434 0.868
## 10 33Hz dispersions 7.870 2 0.020 0.040
## 11 40Hz polar vectors 0.897 2 0.639 1.000
## 12 40Hz dispersions 7.848 2 0.020 0.040
## 13 50Hz polar vectors 0.106 2 0.948 1.000
## 14 50Hz dispersions 2.521 2 0.284 0.568
## 15 80Hz polar vectors 0.806 2 0.668 1.000
## 16 80Hz dispersions 0.492 2 0.782 1.000
rao_testsc <- do.call(rbind, lapply(sort(unique(circ_data_df$condition)), function(x) {
df <- circ_data_df[circ_data_df$condition==x,] %>% group_split(modulation)
rao <- circular::rao.test(lapply(df,'[[', 'phase'))
data.frame(condition = as.factor(x),
test = as.factor(c("polar vectors", "dispersions")),
statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>%
mutate(across(where(is.numeric), round, 3)) %>%
mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value")
}))
# test differences for mod rates within conditions
rao_testsc
## condition test statistic df p_value p_adjust
## 1 silence polar vectors 1.220 7 0.990 1.000
## 2 silence dispersions 1.016 7 0.995 1.000
## 3 steady-state masker polar vectors 36.305 7 0.000 0.000
## 4 steady-state masker dispersions 233.993 7 0.000 0.000
## 5 random masker polar vectors 13.431 7 0.062 0.124
## 6 random masker dispersions 277.880 7 0.000 0.000
Next, let’s compute post-hoc comparisons for each significant test (by modulation rate) above by comparing each pair of masking conditions within each modulation condition for both rates, again using the Rao test.
do_contrasts <- rao_testsm$modulation[rao_testsm$p_value<0.05]
circ_contrasts <- lapply(do_contrasts, function(x) {
df <- circ_data_df[circ_data_df$modulation==x,]
pairs <- t(combn(unique(df$condition), m=2))
tests <- lapply(seq.int(nrow(pairs)), function(y) {
condlist <- df[df$condition %in% pairs[y,],] %>% group_split(condition)
rao <- circular::rao.test(condlist[[1]]$phase, condlist[[2]]$phase)
data.frame(modulation = as.factor(x),
condition = as.factor(paste(pairs[y,],collapse = " vs. ")),
test = as.factor(c("polar vectors", "dispersions")),
statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>%
mutate(across(is.numeric, round, 3)) %>%
mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value")
}) %>% do.call(rbind,.)
}
) %>% do.call(rbind,.)
circ_contrasts
## modulation condition test statistic df p_value p_adjust
## 1 4Hz silence vs. steady-state masker polar vectors 2.985 1 0.084 0.168
## 2 4Hz silence vs. steady-state masker dispersions 79.850 1 0.000 0.000
## 3 4Hz silence vs. random masker polar vectors 3.087 1 0.079 0.158
## 4 4Hz silence vs. random masker dispersions 153.051 1 0.000 0.000
## 5 4Hz steady-state masker vs. random masker polar vectors 0.000 1 0.983 1.000
## 6 4Hz steady-state masker vs. random masker dispersions 15.231 1 0.000 0.000
## 7 8Hz silence vs. steady-state masker polar vectors 0.046 1 0.830 1.000
## 8 8Hz silence vs. steady-state masker dispersions 104.750 1 0.000 0.000
## 9 8Hz silence vs. random masker polar vectors 0.000 1 0.985 1.000
## 10 8Hz silence vs. random masker dispersions 95.520 1 0.000 0.000
## 11 8Hz steady-state masker vs. random masker polar vectors 18.858 1 0.000 0.000
## 12 8Hz steady-state masker vs. random masker dispersions 0.944 1 0.331 0.662
## 13 8Hz silence vs. steady-state masker polar vectors 0.046 1 0.830 1.000
## 14 8Hz silence vs. steady-state masker dispersions 104.750 1 0.000 0.000
## 15 8Hz silence vs. random masker polar vectors 0.000 1 0.985 1.000
## 16 8Hz silence vs. random masker dispersions 95.520 1 0.000 0.000
## 17 8Hz steady-state masker vs. random masker polar vectors 18.858 1 0.000 0.000
## 18 8Hz steady-state masker vs. random masker dispersions 0.944 1 0.331 0.662
## 19 16Hz silence vs. steady-state masker polar vectors 0.848 1 0.357 0.714
## 20 16Hz silence vs. steady-state masker dispersions 51.665 1 0.000 0.000
## 21 16Hz silence vs. random masker polar vectors 1.579 1 0.209 0.418
## 22 16Hz silence vs. random masker dispersions 25.990 1 0.000 0.000
## 23 16Hz steady-state masker vs. random masker polar vectors 5.969 1 0.015 0.030
## 24 16Hz steady-state masker vs. random masker dispersions 0.470 1 0.493 0.986
## 25 16Hz silence vs. steady-state masker polar vectors 0.848 1 0.357 0.714
## 26 16Hz silence vs. steady-state masker dispersions 51.665 1 0.000 0.000
## 27 16Hz silence vs. random masker polar vectors 1.579 1 0.209 0.418
## 28 16Hz silence vs. random masker dispersions 25.990 1 0.000 0.000
## 29 16Hz steady-state masker vs. random masker polar vectors 5.969 1 0.015 0.030
## 30 16Hz steady-state masker vs. random masker dispersions 0.470 1 0.493 0.986
## 31 25Hz silence vs. steady-state masker polar vectors 0.044 1 0.834 1.000
## 32 25Hz silence vs. steady-state masker dispersions 16.873 1 0.000 0.000
## 33 25Hz silence vs. random masker polar vectors 0.008 1 0.927 1.000
## 34 25Hz silence vs. random masker dispersions 5.805 1 0.016 0.032
## 35 25Hz steady-state masker vs. random masker polar vectors 3.840 1 0.050 0.100
## 36 25Hz steady-state masker vs. random masker dispersions 0.732 1 0.392 0.784
## 37 33Hz silence vs. steady-state masker polar vectors 1.651 1 0.199 0.398
## 38 33Hz silence vs. steady-state masker dispersions 0.629 1 0.428 0.856
## 39 33Hz silence vs. random masker polar vectors 0.016 1 0.899 1.000
## 40 33Hz silence vs. random masker dispersions 7.305 1 0.007 0.014
## 41 33Hz steady-state masker vs. random masker polar vectors 0.018 1 0.892 1.000
## 42 33Hz steady-state masker vs. random masker dispersions 6.521 1 0.011 0.022
## 43 40Hz silence vs. steady-state masker polar vectors 0.849 1 0.357 0.714
## 44 40Hz silence vs. steady-state masker dispersions 6.625 1 0.010 0.020
## 45 40Hz silence vs. random masker polar vectors 0.304 1 0.582 1.000
## 46 40Hz silence vs. random masker dispersions 1.231 1 0.267 0.534
## 47 40Hz steady-state masker vs. random masker polar vectors 0.289 1 0.591 1.000
## 48 40Hz steady-state masker vs. random masker dispersions 0.362 1 0.547 1.000
Let’s look at the overall number of calls observed in the three experimental conditions and within the two tested modulation rates.
Two post-hoc hypotheses: 1. The number of calls observed in the silent condition would be be greater than that in the masking conditions. - Method: Pearson’s Chi-squared Test for Goodness-of-Fit - H1: Observed calls unevenly distributed across groups.
Contingency tables:
n_calls <- data %>% group_by(modulation, condition) %>% summarise(n = n())
(n_conttable <- addmargins(table(data$modulation,data$condition), c(1,2)))
##
## silence steady-state masker random masker Sum
## 4Hz 59635 28938 105990 194563
## 8Hz 62203 65446 51939 179588
## 16Hz 65769 60482 25543 151794
## 25Hz 67403 32743 16600 116746
## 33Hz 78179 34843 12830 125852
## 40Hz 68797 34176 10531 113504
## 50Hz 61956 25425 8487 95868
## 80Hz 75144 30675 5260 111079
## Sum 539086 312728 237180 1088994
Goodness of fit tests (test of imbalanced observations) for conditions within each modulation rate:
## 4Hz : X(2, N = 194563) = 46401.95, p < 0.01 (2-tailed)
## 8Hz : X(2, N = 179588) = 1661.06, p < 0.01 (2-tailed)
## 16Hz : X(2, N = 151794) = 18886.24, p < 0.01 (2-tailed)
## 25Hz : X(2, N = 116746) = 34629.51, p < 0.01 (2-tailed)
## 33Hz : X(2, N = 125852) = 52705.34, p < 0.01 (2-tailed)
## 40Hz : X(2, N = 113504) = 45395.99, p < 0.01 (2-tailed)
## 50Hz : X(2, N = 95868) = 46734.51, p < 0.01 (2-tailed)
## 80Hz : X(2, N = 111079) = 67584.22, p < 0.01 (2-tailed)
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
Goodness of fit tests for rates within conditions:
## silence : X(7, N = 539086) = 4417.96, p < 0.01 (2-tailed)
## steady-state masker : X(7, N = 312728) = 40810.74, p < 0.01 (2-tailed)
## random masker : X(7, N = 237180) = 276682.97, p < 0.01 (2-tailed)
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
In the chi-square results below, columns refer to the following data:
observed: the observed counts,expected: the expected counts under the null
hypothesis,residuals: the Pearson residuals, \((observed - expected) /
\sqrt(expected)\),stdres: standardized residuals, \((observed - expected) / \sqrt(V)\),where \(V\) is the residual cell variance.
## Pearson's Chi-squared test for each modulation rate x condition:
## Mods vs. Conds: X(14, N = 1088994) = 209611.08, p < 0.01 (2-tailed)
## # A tibble: 24 × 9
## condition modulation observed prop row_prop col_prop expected resid std_resid
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 silence 4Hz 59635 0.055 0.111 0.307 96315. -118. -184.
## 2 steady-state masker 4Hz 28938 0.027 0.093 0.149 55873. -114. -149.
## 3 random masker 4Hz 105990 0.097 0.447 0.545 42375. 309. 386.
## 4 silence 8Hz 62203 0.057 0.115 0.346 88902. -89.5 -138.
## 5 steady-state masker 8Hz 65446 0.06 0.209 0.364 51573. 61.1 79.2
## 6 random masker 8Hz 51939 0.048 0.219 0.289 39114. 64.8 80.2
## 7 silence 16Hz 65769 0.06 0.122 0.433 75143. -34.2 -51.9
## 8 steady-state masker 16Hz 60482 0.056 0.193 0.398 43591. 80.9 103.
## 9 random masker 16Hz 25543 0.023 0.108 0.168 33060. -41.3 -50.4
## 10 silence 25Hz 67403 0.062 0.125 0.577 57793. 40.0 59.5
## # … with 14 more rows
First we modelled our counts using a Poisson GLM. However, we observed that the data were highly overdispersed. For example, for the 8 Hz context, overdispersion was:
# Poisson Regression
calls.glm <- glm(n ~ condition, family = poisson, data = n_calls_reg[n_calls_reg$modulation=="8Hz",])
# model summary
# summary(calls.glm)
# check for overdispersion
OD.pr <- sum(residuals(calls.glm, type="pearson")^2)/df.residual(calls.glm)
print(OD.pr)
## [1] 6718.167
Therefore, we proceeded with a Negative Binomial GLM, which allows for overdispersion in the count variable.
We ran a model for each modulation rate separately, so that the reference group (Intercept) is meaningful for all model terms for our purposes, i.e. coefficient indicate changes in calling rate between silent baseline and each masking condition.
Our models are given by: \[ln(\widehat{calls_{i}}) = Intercept + \beta_1 Mask(condition_j = 2) + \beta_3 Mask(condition_j = 3)\]
Where \(i\) gives the modulation rate and \(j\) gives the level of the condition manipulation.
Thus, \[\widehat{calls_i} = e^{Intercept + \beta_1 Mask_{condition_j = 2}+...}\] defines the rate of calling per hour (the sampling time for each condition), given an arbitrary combination of modulation rate and masking condition.
We also ran a “big” model that included all interactions, for which coefficients would have a more complex meaning but analyses of variance would provide a global picture.
big.model <- run_negbinom_model(n_calls_reg,
formula = "n ~ condition * modulation",
contformula = "~ modulation | condition ")
## [1] "~~~~~~~~~~~~~~ all rates ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.4410757978,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4987 -1.4098 -0.5586 0.2112 1.8856
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.00027 0.33671 23.760 < 0.0000000000000002 ***
## conditionsteady-state masker -0.72309 0.47620 -1.518 0.128903
## conditionrandom masker 0.57510 0.47618 1.208 0.227144
## modulation8Hz 0.04216 0.47618 0.089 0.929449
## modulation16Hz 0.09791 0.47618 0.206 0.837099
## modulation25Hz 0.12245 0.47618 0.257 0.797069
## modulation33Hz 0.27076 0.47618 0.569 0.569625
## modulation40Hz 0.19421 0.48241 0.403 0.687252
## modulation50Hz 0.03818 0.47618 0.080 0.936092
## modulation80Hz 0.23116 0.47618 0.485 0.627355
## conditionsteady-state masker:modulation8Hz 0.82520 0.67785 1.217 0.223461
## conditionrandom masker:modulation8Hz -0.75544 0.67342 -1.122 0.261954
## conditionsteady-state masker:modulation16Hz 0.63928 0.67344 0.949 0.342477
## conditionrandom masker:modulation16Hz -1.52089 0.67344 -2.258 0.023921 *
## conditionsteady-state masker:modulation25Hz 0.05238 0.67786 0.077 0.938407
## conditionrandom masker:modulation25Hz -1.97639 0.67345 -2.935 0.003339 **
## conditionsteady-state masker:modulation33Hz -0.08506 0.67345 -0.126 0.899488
## conditionrandom masker:modulation33Hz -2.38232 0.67346 -3.537 0.000404 ***
## conditionsteady-state masker:modulation40Hz 0.02345 0.68225 0.034 0.972580
## conditionrandom masker:modulation40Hz -2.50323 0.67789 -3.693 0.000222 ***
## conditionsteady-state masker:modulation50Hz -0.16760 0.67346 -0.249 0.803460
## conditionrandom masker:modulation50Hz -2.56299 0.67350 -3.806 0.000142 ***
## conditionsteady-state masker:modulation80Hz -0.17287 0.67345 -0.257 0.797414
## conditionrandom masker:modulation80Hz -3.23438 0.67355 -4.802 0.00000157 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4411) family taken to be 1)
##
## Null deviance: 713.90 on 475 degrees of freedom
## Residual deviance: 613.84 on 452 degrees of freedom
## AIC: 7837
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4411
## Std. Err.: 0.0234
##
## 2 x log-likelihood: -7786.9700
## [1] "Model dispersion..."
## [1] 0.85917
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 45.676 2 0.0000000001207 ***
## modulation 32.121 7 0.0000385773970 ***
## condition:modulation 42.692 14 0.0000959364219 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 2981.75000000 1645.30186853 6257.4686715
## conditionsteady-state masker 0.48525195 0.18758610 1.2552607
## conditionrandom masker 1.77731198 0.68709741 4.5973654
## modulation8Hz 1.04306196 0.40323584 2.6981189
## modulation16Hz 1.10285906 0.42635344 2.8527930
## modulation25Hz 1.13025908 0.43694631 2.9236672
## modulation33Hz 1.31095833 0.50680476 3.3910726
## modulation40Hz 1.21435222 0.46526044 3.1943692
## modulation50Hz 1.03892010 0.40163459 2.6874054
## modulation80Hz 1.26006540 0.48712952 3.2594305
## conditionsteady-state masker:modulation8Hz 2.28234285 0.59826840 8.7411771
## conditionrandom masker:modulation8Hz 0.46980603 0.12401175 1.7798129
## conditionsteady-state masker:modulation16Hz 1.89512393 0.50022870 7.1797054
## conditionrandom masker:modulation16Hz 0.21851789 0.05767924 0.8278554
## conditionsteady-state masker:modulation25Hz 1.05377611 0.27621988 4.0359513
## conditionrandom masker:modulation25Hz 0.13856872 0.03657505 0.5249832
## conditionsteady-state masker:modulation33Hz 0.91845555 0.24242820 3.4796307
## conditionrandom masker:modulation33Hz 0.09233639 0.02437154 0.3498346
## conditionsteady-state masker:modulation40Hz 1.02372767 0.26545814 3.9479609
## conditionrandom masker:modulation40Hz 0.08182011 0.02136184 0.3121601
## conditionsteady-state masker:modulation50Hz 0.84568826 0.22321671 3.2040103
## conditionrandom masker:modulation50Hz 0.07707387 0.02034186 0.2920274
## conditionsteady-state masker:modulation80Hz 0.84124592 0.22204717 3.1871367
## conditionrandom masker:modulation80Hz 0.03938472 0.01039366 0.1492406
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition = silence:
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 2982 1004.0 Inf 1541 5769
## 8Hz 3110 1047.2 Inf 1608 6017
## 16Hz 3288 1107.3 Inf 1700 6362
## 25Hz 3370 1134.8 Inf 1742 6520
## 33Hz 3909 1316.2 Inf 2020 7563
## 40Hz 3621 1250.9 Inf 1840 7126
## 50Hz 3098 1043.1 Inf 1601 5993
## 80Hz 3757 1265.1 Inf 1942 7269
##
## condition = steady-state masker:
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 1447 487.2 Inf 748 2799
## 8Hz 3445 1189.9 Inf 1750 6779
## 16Hz 3024 1018.3 Inf 1563 5851
## 25Hz 1723 595.4 Inf 876 3392
## 33Hz 1742 586.6 Inf 900 3371
## 40Hz 1799 621.4 Inf 914 3540
## 50Hz 1271 428.1 Inf 657 2460
## 80Hz 1534 516.5 Inf 793 2967
##
## condition = random masker:
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 5300 1784.4 Inf 2739 10253
## 8Hz 2597 874.4 Inf 1342 5024
## 16Hz 1277 430.1 Inf 660 2471
## 25Hz 830 279.5 Inf 429 1606
## 33Hz 642 216.1 Inf 332 1241
## 40Hz 527 177.4 Inf 272 1019
## 50Hz 424 142.9 Inf 219 821
## 80Hz 263 88.6 Inf 136 509
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## condition = silence:
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 0.959 0.457 Inf 1 -0.089 1.0000
## 4Hz / 16Hz 0.907 0.432 Inf 1 -0.206 1.0000
## 4Hz / 25Hz 0.885 0.421 Inf 1 -0.257 1.0000
## 4Hz / 33Hz 0.763 0.363 Inf 1 -0.569 1.0000
## 4Hz / 40Hz 0.823 0.397 Inf 1 -0.403 1.0000
## 4Hz / 50Hz 0.963 0.458 Inf 1 -0.080 1.0000
## 4Hz / 80Hz 0.794 0.378 Inf 1 -0.485 1.0000
## 8Hz / 16Hz 0.946 0.450 Inf 1 -0.117 1.0000
## 8Hz / 25Hz 0.923 0.439 Inf 1 -0.169 1.0000
## 8Hz / 33Hz 0.796 0.379 Inf 1 -0.480 1.0000
## 8Hz / 40Hz 0.859 0.414 Inf 1 -0.315 1.0000
## 8Hz / 50Hz 1.004 0.478 Inf 1 0.008 1.0000
## 8Hz / 80Hz 0.828 0.394 Inf 1 -0.397 1.0000
## 16Hz / 25Hz 0.976 0.465 Inf 1 -0.052 1.0000
## 16Hz / 33Hz 0.841 0.401 Inf 1 -0.363 1.0000
## 16Hz / 40Hz 0.908 0.438 Inf 1 -0.200 1.0000
## 16Hz / 50Hz 1.062 0.505 Inf 1 0.125 1.0000
## 16Hz / 80Hz 0.875 0.417 Inf 1 -0.280 1.0000
## 25Hz / 33Hz 0.862 0.411 Inf 1 -0.311 1.0000
## 25Hz / 40Hz 0.931 0.449 Inf 1 -0.149 1.0000
## 25Hz / 50Hz 1.088 0.518 Inf 1 0.177 1.0000
## 25Hz / 80Hz 0.897 0.427 Inf 1 -0.228 1.0000
## 33Hz / 40Hz 1.080 0.521 Inf 1 0.159 1.0000
## 33Hz / 50Hz 1.262 0.601 Inf 1 0.488 1.0000
## 33Hz / 80Hz 1.040 0.495 Inf 1 0.083 1.0000
## 40Hz / 50Hz 1.169 0.564 Inf 1 0.323 1.0000
## 40Hz / 80Hz 0.964 0.465 Inf 1 -0.077 1.0000
## 50Hz / 80Hz 0.824 0.393 Inf 1 -0.405 1.0000
##
## condition = steady-state masker:
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 0.420 0.203 Inf 1 -1.798 1.0000
## 4Hz / 16Hz 0.478 0.228 Inf 1 -1.548 1.0000
## 4Hz / 25Hz 0.840 0.405 Inf 1 -0.362 1.0000
## 4Hz / 33Hz 0.831 0.396 Inf 1 -0.390 1.0000
## 4Hz / 40Hz 0.804 0.388 Inf 1 -0.451 1.0000
## 4Hz / 50Hz 1.138 0.542 Inf 1 0.272 1.0000
## 4Hz / 80Hz 0.943 0.449 Inf 1 -0.122 1.0000
## 8Hz / 16Hz 1.139 0.549 Inf 1 0.270 1.0000
## 8Hz / 25Hz 1.999 0.977 Inf 1 1.417 1.0000
## 8Hz / 33Hz 1.977 0.954 Inf 1 1.413 1.0000
## 8Hz / 40Hz 1.915 0.936 Inf 1 1.330 1.0000
## 8Hz / 50Hz 2.710 1.307 Inf 1 2.066 1.0000
## 8Hz / 80Hz 2.246 1.083 Inf 1 1.677 1.0000
## 16Hz / 25Hz 1.755 0.847 Inf 1 1.166 1.0000
## 16Hz / 33Hz 1.736 0.827 Inf 1 1.158 1.0000
## 16Hz / 40Hz 1.681 0.811 Inf 1 1.077 1.0000
## 16Hz / 50Hz 2.379 1.133 Inf 1 1.820 1.0000
## 16Hz / 80Hz 1.972 0.939 Inf 1 1.426 1.0000
## 25Hz / 33Hz 0.989 0.477 Inf 1 -0.023 1.0000
## 25Hz / 40Hz 0.958 0.468 Inf 1 -0.088 1.0000
## 25Hz / 50Hz 1.356 0.654 Inf 1 0.631 1.0000
## 25Hz / 80Hz 1.124 0.542 Inf 1 0.242 1.0000
## 33Hz / 40Hz 0.969 0.467 Inf 1 -0.066 1.0000
## 33Hz / 50Hz 1.370 0.653 Inf 1 0.662 1.0000
## 33Hz / 80Hz 1.136 0.541 Inf 1 0.268 1.0000
## 40Hz / 50Hz 1.415 0.683 Inf 1 0.719 1.0000
## 40Hz / 80Hz 1.173 0.566 Inf 1 0.330 1.0000
## 50Hz / 80Hz 0.829 0.395 Inf 1 -0.394 1.0000
##
## condition = random masker:
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 2.041 0.972 Inf 1 1.498 1.0000
## 4Hz / 16Hz 4.149 1.976 Inf 1 2.988 0.0786
## 4Hz / 25Hz 6.385 3.041 Inf 1 3.893 0.0028
## 4Hz / 33Hz 8.261 3.934 Inf 1 4.434 0.0003
## 4Hz / 40Hz 10.065 4.793 Inf 1 4.848 <.0001
## 4Hz / 50Hz 12.489 5.948 Inf 1 5.301 <.0001
## 4Hz / 80Hz 20.150 9.599 Inf 1 6.305 <.0001
## 8Hz / 16Hz 2.033 0.968 Inf 1 1.490 1.0000
## 8Hz / 25Hz 3.129 1.490 Inf 1 2.395 0.4651
## 8Hz / 33Hz 4.048 1.928 Inf 1 2.936 0.0931
## 8Hz / 40Hz 4.932 2.349 Inf 1 3.351 0.0226
## 8Hz / 50Hz 6.120 2.915 Inf 1 3.803 0.0040
## 8Hz / 80Hz 9.874 4.704 Inf 1 4.807 <.0001
## 16Hz / 25Hz 1.539 0.733 Inf 1 0.905 1.0000
## 16Hz / 33Hz 1.991 0.948 Inf 1 1.446 1.0000
## 16Hz / 40Hz 2.426 1.155 Inf 1 1.860 1.0000
## 16Hz / 50Hz 3.010 1.434 Inf 1 2.313 0.5799
## 16Hz / 80Hz 4.856 2.313 Inf 1 3.317 0.0255
## 25Hz / 33Hz 1.294 0.616 Inf 1 0.541 1.0000
## 25Hz / 40Hz 1.576 0.751 Inf 1 0.955 1.0000
## 25Hz / 50Hz 1.956 0.932 Inf 1 1.408 1.0000
## 25Hz / 80Hz 3.156 1.504 Inf 1 2.412 0.4438
## 33Hz / 40Hz 1.218 0.580 Inf 1 0.415 1.0000
## 33Hz / 50Hz 1.512 0.720 Inf 1 0.868 1.0000
## 33Hz / 80Hz 2.439 1.162 Inf 1 1.872 1.0000
## 40Hz / 50Hz 1.241 0.591 Inf 1 0.453 1.0000
## 40Hz / 80Hz 2.002 0.954 Inf 1 1.457 1.0000
## 50Hz / 80Hz 1.613 0.769 Inf 1 1.004 1.0000
##
## P value adjustment: bonferroni method for 28 tests
## Tests are performed on the log scale
# overdispersion
print(big.model$dispersion)
## [1] 0.85917
all.models <- lapply(sort(unique(n_calls_reg$modulation)), function(x)
md <- run_negbinom_model(n_calls_reg %>% filter(modulation==x)) )
## [1] "~~~~~~~~~~~~~~ 4Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4723339773,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3464 -1.3494 -0.7116 0.1990 1.5790
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0003 0.3254 24.587 <0.0000000000000002 ***
## conditionsteady-state masker -0.7231 0.4602 -1.571 0.116
## conditionrandom masker 0.5751 0.4602 1.250 0.211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4723) family taken to be 1)
##
## Null deviance: 84.212 on 59 degrees of freedom
## Residual deviance: 76.661 on 57 degrees of freedom
## AIC: 1051.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4723
## Std. Err.: 0.0710
##
## 2 x log-likelihood: -1043.1210
## [1] "Model dispersion..."
## [1] 0.9538492
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 7.5502 2 0.02294 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 2981.7500000 1675.1885512 6085.738182
## conditionsteady-state masker 0.4852519 0.1938973 1.214403
## conditionrandom masker 1.7773120 0.7102157 4.447716
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 2982 970 Inf 1576 5642
## steady-state masker 1447 471 Inf 765 2738
## random masker 5300 1724 Inf 2801 10027
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 2.061 0.948 Inf 1 1.571 0.3483
## silence / random masker 0.563 0.259 Inf 1 -1.250 0.6341
## (steady-state masker) / random masker 0.273 0.126 Inf 1 -2.821 0.0144
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 8Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.3868466812,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.32263 -1.40195 -0.71589 0.06566 1.42301
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0424 0.3595 22.369 <0.0000000000000002 ***
## conditionsteady-state masker 0.1021 0.5151 0.198 0.843
## conditionrandom masker -0.1803 0.5085 -0.355 0.723
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3868) family taken to be 1)
##
## Null deviance: 77.915 on 58 degrees of freedom
## Residual deviance: 77.609 on 56 degrees of freedom
## AIC: 1017.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3868
## Std. Err.: 0.0577
##
## 2 x log-likelihood: -1009.7670
## [1] "Model dispersion..."
## [1] 0.8545476
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 0.3059 2 0.8582
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3110.1500000 1655.2222327 6909.314899
## conditionsteady-state masker 1.1075113 0.3966732 3.119991
## conditionrandom masker 0.8349919 0.3018830 2.309542
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3110 1118 Inf 1537 6292
## steady-state masker 3445 1271 Inf 1672 7098
## random masker 2597 934 Inf 1284 5254
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 0.903 0.465 Inf 1 -0.198 1.0000
## silence / random masker 1.198 0.609 Inf 1 0.355 1.0000
## (steady-state masker) / random masker 1.326 0.683 Inf 1 0.548 1.0000
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 16Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.5221485665,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2929 -1.2910 -0.4686 0.3407 1.5296
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0982 0.3095 26.168 <0.0000000000000002 ***
## conditionsteady-state masker -0.0838 0.4377 -0.191 0.8482
## conditionrandom masker -0.9458 0.4377 -2.161 0.0307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.5221) family taken to be 1)
##
## Null deviance: 80.629 on 59 degrees of freedom
## Residual deviance: 75.558 on 57 degrees of freedom
## AIC: 1036.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.5221
## Std. Err.: 0.0792
##
## 2 x log-likelihood: -1028.4460
## [1] "Model dispersion..."
## [1] 0.8309316
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 5.0715 2 0.0792 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3288.4500000 1895.4786790 6456.3843402
## conditionsteady-state masker 0.9196126 0.3848699 2.1973327
## conditionrandom masker 0.3883745 0.1625319 0.9280312
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3288 1018 Inf 1793 6031
## steady-state masker 3024 936 Inf 1649 5547
## random masker 1277 395 Inf 696 2343
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 1.09 0.476 Inf 1 0.191 1.0000
## silence / random masker 2.57 1.127 Inf 1 2.161 0.0921
## (steady-state masker) / random masker 2.37 1.036 Inf 1 1.969 0.1467
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 25Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4039458612,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2513 -1.3790 -0.7238 0.1119 1.3833
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.1227 0.3518 23.086 < 0.0000000000000002 ***
## conditionsteady-state masker -0.6707 0.5041 -1.331 0.18335
## conditionrandom masker -1.4013 0.4976 -2.816 0.00486 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4039) family taken to be 1)
##
## Null deviance: 84.625 on 58 degrees of freedom
## Residual deviance: 77.072 on 56 degrees of freedom
## AIC: 954.96
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4039
## Std. Err.: 0.0604
##
## 2 x log-likelihood: -946.9580
## [1] "Model dispersion..."
## [1] 0.8618375
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 7.5536 2 0.0229 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3370.1500000 1815.42972059 7344.0187496
## conditionsteady-state masker 0.5113469 0.18736028 1.4075862
## conditionrandom masker 0.2462798 0.09107138 0.6660024
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3370 1186 Inf 1691 6716
## steady-state masker 1723 622 Inf 849 3497
## random masker 830 292 Inf 416 1654
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 1.96 0.986 Inf 1 1.331 0.5501
## silence / random masker 4.06 2.021 Inf 1 2.816 0.0146
## (steady-state masker) / random masker 2.08 1.047 Inf 1 1.449 0.4419
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 33Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4259724824,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4566 -1.3961 -0.5122 0.2072 1.8530
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.2710 0.3426 24.140 < 0.0000000000000002 ***
## conditionsteady-state masker -0.8081 0.4846 -1.668 0.095356 .
## conditionrandom masker -1.8072 0.4846 -3.729 0.000192 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.426) family taken to be 1)
##
## Null deviance: 90.429 on 59 degrees of freedom
## Residual deviance: 77.690 on 57 degrees of freedom
## AIC: 973.91
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4260
## Std. Err.: 0.0635
##
## 2 x log-likelihood: -965.9070
## [1] "Model dispersion..."
## [1] 0.8952745
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 12.739 2 0.001713 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3908.9500000 2136.57666836 8324.8275436
## conditionsteady-state masker 0.4456823 0.16933647 1.1730063
## conditionrandom masker 0.1641106 0.06234774 0.4319688
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3909 1339 Inf 1997 7651
## steady-state masker 1742 597 Inf 890 3410
## random masker 642 220 Inf 328 1256
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 2.24 1.09 Inf 1 1.668 0.2861
## silence / random masker 6.09 2.95 Inf 1 3.729 0.0006
## (steady-state masker) / random masker 2.72 1.32 Inf 1 2.062 0.1178
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 40Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4689172885,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0805 -1.4607 -0.5346 0.3031 1.4668
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.1945 0.3350 24.458 < 0.0000000000000002 ***
## conditionsteady-state masker -0.6996 0.4738 -1.477 0.14
## conditionrandom masker -1.9281 0.4680 -4.120 0.0000378 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4689) family taken to be 1)
##
## Null deviance: 89.353 on 57 degrees of freedom
## Residual deviance: 74.197 on 55 degrees of freedom
## AIC: 940.19
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4689
## Std. Err.: 0.0716
##
## 2 x log-likelihood: -932.1910
## [1] "Model dispersion..."
## [1] 0.841219
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 15.156 2 0.0005115 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3620.8947368 2003.10860148 7568.0222469
## conditionsteady-state masker 0.4967658 0.19297987 1.2787671
## conditionrandom masker 0.1454199 0.05697374 0.3684498
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3621 1213 Inf 1878 6982
## steady-state masker 1799 603 Inf 933 3469
## random masker 527 172 Inf 278 999
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 2.01 0.954 Inf 1 1.477 0.4194
## silence / random masker 6.88 3.218 Inf 1 4.120 0.0001
## (steady-state masker) / random masker 3.42 1.599 Inf 1 2.625 0.0260
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 50Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4025397372,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.17513 -1.47676 -0.69721 0.04592 1.42237
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0384 0.3525 22.807 < 0.0000000000000002 ***
## conditionsteady-state masker -0.8907 0.4985 -1.787 0.074 .
## conditionrandom masker -1.9879 0.4986 -3.987 0.0000668 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4025) family taken to be 1)
##
## Null deviance: 92.781 on 59 degrees of freedom
## Residual deviance: 78.417 on 57 degrees of freedom
## AIC: 928.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4025
## Std. Err.: 0.0597
##
## 2 x log-likelihood: -920.2040
## [1] "Model dispersion..."
## [1] 0.8204445
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 14.364 2 0.0007603 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3097.8000000 1667.09359936 6760.9100740
## conditionsteady-state masker 0.4103719 0.15148246 1.1117133
## conditionrandom masker 0.1369843 0.05055829 0.3711498
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3098 1092 Inf 1553 6181
## steady-state masker 1271 448 Inf 637 2537
## random masker 424 150 Inf 213 847
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 2.44 1.21 Inf 1 1.787 0.2219
## silence / random masker 7.30 3.64 Inf 1 3.987 0.0002
## (steady-state masker) / random masker 3.00 1.49 Inf 1 2.201 0.0833
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ 80Hz ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4827242071,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0092 -1.3175 -0.5869 0.2676 1.4456
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.2314 0.3219 25.575 < 0.0000000000000002 ***
## conditionsteady-state masker -0.8960 0.4552 -1.968 0.049 *
## conditionrandom masker -2.6593 0.4554 -5.840 0.00000000523 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4827) family taken to be 1)
##
## Null deviance: 103.877 on 59 degrees of freedom
## Residual deviance: 76.242 on 57 degrees of freedom
## AIC: 944.88
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4827
## Std. Err.: 0.0728
##
## 2 x log-likelihood: -936.8780
## [1] "Model dispersion..."
## [1] 0.8198301
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## condition 27.634 2 0.0000009984 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 3757.20000000 2122.8400904 7602.5475769
## conditionsteady-state masker 0.40821622 0.1647983 1.0111782
## conditionrandom masker 0.06999894 0.0282497 0.1734479
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition response SE df asymp.LCL asymp.UCL
## silence 3757 1209.3 Inf 1999 7060
## steady-state masker 1534 493.7 Inf 816 2882
## random masker 263 84.7 Inf 140 494
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## silence / (steady-state masker) 2.45 1.12 Inf 1 1.968 0.1471
## silence / random masker 14.29 6.51 Inf 1 5.840 <.0001
## (steady-state masker) / random masker 5.83 2.66 Inf 1 3.872 0.0003
##
## P value adjustment: bonferroni method for 3 tests
## Tests are performed on the log scale
names(all.models) <- sort(unique(n_calls_reg$modulation))
# overdispersion
print(do.call(rbind, lapply(all.models, '[[', 'dispersion')) )
## [,1]
## 4Hz 0.9538492
## 8Hz 0.8545476
## 16Hz 0.8309316
## 25Hz 0.8618375
## 33Hz 0.8952745
## 40Hz 0.8412190
## 50Hz 0.8204445
## 80Hz 0.8198301
Note that these are 1-way ANOVAs from different models put together in one table, not a single ANOVA.
# pull out anova tables
(anova <- do.call(rbind, lapply(all.models, '[[', 'deviance')) %>%
mutate(across(is.numeric, round, 2)))
## LR Chisq Df Pr(>Chisq)
## 4Hz 7.55 2 0.02
## 8Hz 0.31 2 0.86
## 16Hz 5.07 2 0.08
## 25Hz 7.55 2 0.02
## 33Hz 12.74 2 0.00
## 40Hz 15.16 2 0.00
## 50Hz 14.36 2 0.00
## 80Hz 27.63 2 0.00
Below, the analysis of variance for the large model with all features:
big.model$deviance
## LR Chisq Df Pr(>Chisq)
## condition 45.676 2 0
## modulation 32.121 7 0
## condition:modulation 42.692 14 0
The IRRs give us the change in the per hour rate [of calling] from baseline to each condition.
# IRR
IRRs <- do.call(rbind, lapply(all.models, '[[', 'IRR')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(modulation = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
mutate(coef = str_remove(str_split(row.names(.),"[.]",simplify = T)[,2],"condition"), .after=1) %>%
remove_rownames()
IRRs
## modulation coef Estimate 2.5 % 97.5 %
## 1 4Hz (Intercept) 2981.75 1675.19 6085.74
## 2 4Hz steady-state masker 0.49 0.19 1.21
## 3 4Hz random masker 1.78 0.71 4.45
## 4 8Hz (Intercept) 3110.15 1655.22 6909.31
## 5 8Hz steady-state masker 1.11 0.40 3.12
## 6 8Hz random masker 0.83 0.30 2.31
## 7 16Hz (Intercept) 3288.45 1895.48 6456.38
## 8 16Hz steady-state masker 0.92 0.38 2.20
## 9 16Hz random masker 0.39 0.16 0.93
## 10 25Hz (Intercept) 3370.15 1815.43 7344.02
## 11 25Hz steady-state masker 0.51 0.19 1.41
## 12 25Hz random masker 0.25 0.09 0.67
## 13 33Hz (Intercept) 3908.95 2136.58 8324.83
## 14 33Hz steady-state masker 0.45 0.17 1.17
## 15 33Hz random masker 0.16 0.06 0.43
## 16 40Hz (Intercept) 3620.89 2003.11 7568.02
## 17 40Hz steady-state masker 0.50 0.19 1.28
## 18 40Hz random masker 0.15 0.06 0.37
## 19 50Hz (Intercept) 3097.80 1667.09 6760.91
## 20 50Hz steady-state masker 0.41 0.15 1.11
## 21 50Hz random masker 0.14 0.05 0.37
## 22 80Hz (Intercept) 3757.20 2122.84 7602.55
## 23 80Hz steady-state masker 0.41 0.16 1.01
## 24 80Hz random masker 0.07 0.03 0.17
Below are the post-hoc contrasts from estimated marginal means:
(contrasts <- do.call(rbind, lapply(all.models, '[[', 'contrasts')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(modulation = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
remove_rownames() %>%
dplyr::select(-c("df", "null")) %>%
dplyr::rename(IRR = ratio))
## modulation masking_condition IRR se asymp_lcl asymp_ucl z_ratio p_value
## 1 4Hz silence / (steady-state masker) 2.06 0.95 0.69 6.20 1.57 0.35
## 2 4Hz silence / random masker 0.56 0.26 0.19 1.69 -1.25 0.63
## 3 4Hz (steady-state masker) / random masker 0.27 0.13 0.09 0.82 -2.82 0.01
## 4 8Hz silence / (steady-state masker) 0.90 0.47 0.26 3.10 -0.20 1.00
## 5 8Hz silence / random masker 1.20 0.61 0.36 4.04 0.36 1.00
## 6 8Hz (steady-state masker) / random masker 1.33 0.68 0.39 4.55 0.55 1.00
## 7 16Hz silence / (steady-state masker) 1.09 0.48 0.38 3.10 0.19 1.00
## 8 16Hz silence / random masker 2.58 1.13 0.90 7.34 2.16 0.09
## 9 16Hz (steady-state masker) / random masker 2.37 1.04 0.83 6.75 1.97 0.15
## 10 25Hz silence / (steady-state masker) 1.96 0.99 0.58 6.54 1.33 0.55
## 11 25Hz silence / random masker 4.06 2.02 1.23 13.36 2.82 0.01
## 12 25Hz (steady-state masker) / random masker 2.08 1.05 0.62 6.94 1.45 0.44
## 13 33Hz silence / (steady-state masker) 2.24 1.09 0.70 7.16 1.67 0.29
## 14 33Hz silence / random masker 6.09 2.95 1.91 19.44 3.73 0.00
## 15 33Hz (steady-state masker) / random masker 2.72 1.32 0.85 8.66 2.06 0.12
## 16 40Hz silence / (steady-state masker) 2.01 0.95 0.65 6.26 1.48 0.42
## 17 40Hz silence / random masker 6.88 3.22 2.24 21.08 4.12 0.00
## 18 40Hz (steady-state masker) / random masker 3.42 1.60 1.11 10.47 2.62 0.03
## 19 50Hz silence / (steady-state masker) 2.44 1.22 0.74 8.04 1.79 0.22
## 20 50Hz silence / random masker 7.30 3.64 2.21 24.08 3.99 0.00
## 21 50Hz (steady-state masker) / random masker 3.00 1.49 0.91 9.88 2.20 0.08
## 22 80Hz silence / (steady-state masker) 2.45 1.11 0.82 7.28 1.97 0.15
## 23 80Hz silence / random masker 14.29 6.50 4.80 42.50 5.84 0.00
## 24 80Hz (steady-state masker) / random masker 5.83 2.66 1.96 17.35 3.87 0.00
Now let’s plot the model predictions to ensure it’s a good fit:
preds <- do.call(rbind, lapply(all.models, '[[', 'predictions')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(modulation = factor(str_split(row.names(.),"[.]",simplify = T)[,1],
levels = sort(unique(data$modulation))),.before=1) %>%
remove_rownames()
n_pred <- p.n_pred(preds) %>% recolor()
n_pred
ggsave(file.path(figures_path,paste0(exp, '_fig1e.png') ), n_pred, width = 5, height = 4)
Now let’s see if there were changes in the overall rate of calling between different modulation rate contexts for the same playback condition.
all.models2 <- lapply(sort(unique(n_calls_reg$condition)), function(x)
md <- run_negbinom_model(n_calls_reg %>% filter(condition==x),
formula = "n ~ modulation",
contformula = "~ modulation") )
## [1] "~~~~~~~~~~~~~~ all rates ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.4097400923,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4104 -1.3337 -0.6802 0.2266 1.4645
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.00027 0.34935 22.900 <0.0000000000000002 ***
## modulation8Hz 0.04216 0.49405 0.085 0.932
## modulation16Hz 0.09791 0.49405 0.198 0.843
## modulation25Hz 0.12245 0.49405 0.248 0.804
## modulation33Hz 0.27076 0.49405 0.548 0.584
## modulation40Hz 0.19421 0.50051 0.388 0.698
## modulation50Hz 0.03818 0.49405 0.077 0.938
## modulation80Hz 0.23116 0.49405 0.468 0.640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4097) family taken to be 1)
##
## Null deviance: 208.15 on 158 degrees of freedom
## Residual deviance: 207.59 on 151 degrees of freedom
## AIC: 2794.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4097
## Std. Err.: 0.0374
##
## 2 x log-likelihood: -2776.6580
## [1] "Model dispersion..."
## [1] 0.7875016
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## modulation 0.55954 7 0.9992
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 2981.750000 1612.6030958 6456.973590
## modulation8Hz 1.043062 0.3885751 2.799918
## modulation16Hz 1.102859 0.4108522 2.960428
## modulation25Hz 1.130259 0.4210599 3.033976
## modulation33Hz 1.310958 0.4883784 3.519017
## modulation40Hz 1.214352 0.4482413 3.317702
## modulation50Hz 1.038920 0.3870321 2.788800
## modulation80Hz 1.260065 0.4694185 3.382408
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 2982 1042 Inf 1504 5913
## 8Hz 3110 1087 Inf 1568 6168
## 16Hz 3288 1149 Inf 1658 6522
## 25Hz 3370 1177 Inf 1699 6684
## 33Hz 3909 1366 Inf 1971 7752
## 40Hz 3621 1298 Inf 1794 7310
## 50Hz 3098 1082 Inf 1562 6144
## 80Hz 3757 1313 Inf 1895 7451
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 0.959 0.474 Inf 1 -0.085 1.0000
## 4Hz / 16Hz 0.907 0.448 Inf 1 -0.198 1.0000
## 4Hz / 25Hz 0.885 0.437 Inf 1 -0.248 1.0000
## 4Hz / 33Hz 0.763 0.377 Inf 1 -0.548 1.0000
## 4Hz / 40Hz 0.823 0.412 Inf 1 -0.388 1.0000
## 4Hz / 50Hz 0.963 0.476 Inf 1 -0.077 1.0000
## 4Hz / 80Hz 0.794 0.392 Inf 1 -0.468 1.0000
## 8Hz / 16Hz 0.946 0.467 Inf 1 -0.113 1.0000
## 8Hz / 25Hz 0.923 0.456 Inf 1 -0.163 1.0000
## 8Hz / 33Hz 0.796 0.393 Inf 1 -0.463 1.0000
## 8Hz / 40Hz 0.859 0.430 Inf 1 -0.304 1.0000
## 8Hz / 50Hz 1.004 0.496 Inf 1 0.008 1.0000
## 8Hz / 80Hz 0.828 0.409 Inf 1 -0.383 1.0000
## 16Hz / 25Hz 0.976 0.482 Inf 1 -0.050 1.0000
## 16Hz / 33Hz 0.841 0.416 Inf 1 -0.350 1.0000
## 16Hz / 40Hz 0.908 0.455 Inf 1 -0.192 1.0000
## 16Hz / 50Hz 1.062 0.524 Inf 1 0.121 1.0000
## 16Hz / 80Hz 0.875 0.432 Inf 1 -0.270 1.0000
## 25Hz / 33Hz 0.862 0.426 Inf 1 -0.300 1.0000
## 25Hz / 40Hz 0.931 0.466 Inf 1 -0.143 1.0000
## 25Hz / 50Hz 1.088 0.537 Inf 1 0.171 1.0000
## 25Hz / 80Hz 0.897 0.443 Inf 1 -0.220 1.0000
## 33Hz / 40Hz 1.080 0.540 Inf 1 0.153 1.0000
## 33Hz / 50Hz 1.262 0.623 Inf 1 0.471 1.0000
## 33Hz / 80Hz 1.040 0.514 Inf 1 0.080 1.0000
## 40Hz / 50Hz 1.169 0.585 Inf 1 0.312 1.0000
## 40Hz / 80Hz 0.964 0.482 Inf 1 -0.074 1.0000
## 50Hz / 80Hz 0.824 0.407 Inf 1 -0.391 1.0000
##
## P value adjustment: bonferroni method for 28 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ all rates ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.3986165412,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3569 -1.4586 -0.6550 0.3964 1.7925
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.27718 0.35422 20.545 <0.0000000000000002 ***
## modulation8Hz 0.86736 0.50746 1.709 0.0874 .
## modulation16Hz 0.73719 0.50092 1.472 0.1411
## modulation25Hz 0.17483 0.50748 0.345 0.7305
## modulation33Hz 0.18570 0.50093 0.371 0.7109
## modulation40Hz 0.21766 0.50748 0.429 0.6680
## modulation50Hz -0.12942 0.50094 -0.258 0.7961
## modulation80Hz 0.05829 0.50093 0.116 0.9074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3986) family taken to be 1)
##
## Null deviance: 212.84 on 156 degrees of freedom
## Residual deviance: 205.50 on 149 degrees of freedom
## AIC: 2566.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3986
## Std. Err.: 0.0365
##
## 2 x log-likelihood: -2548.6170
## [1] "Model dispersion..."
## [1] 0.7789167
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## modulation 7.3451 7 0.3939
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 1446.9000000 776.5128310 3171.476821
## modulation8Hz 2.3806250 0.8662196 6.599605
## modulation16Hz 2.0900546 0.7675921 5.690950
## modulation25Hz 1.1910400 0.4333624 3.301910
## modulation33Hz 1.2040569 0.4421913 3.278566
## modulation40Hz 1.2431660 0.4523296 3.446411
## modulation50Hz 0.8786025 0.3226615 2.392422
## modulation80Hz 1.0600249 0.3892925 2.886397
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 1447 513 Inf 723 2897
## 8Hz 3445 1252 Inf 1690 7022
## 16Hz 3024 1071 Inf 1510 6055
## 25Hz 1723 626 Inf 845 3513
## 33Hz 1742 617 Inf 870 3488
## 40Hz 1799 654 Inf 882 3667
## 50Hz 1271 450 Inf 635 2545
## 80Hz 1534 543 Inf 766 3071
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 0.420 0.213 Inf 1 -1.709 1.0000
## 4Hz / 16Hz 0.478 0.240 Inf 1 -1.472 1.0000
## 4Hz / 25Hz 0.840 0.426 Inf 1 -0.345 1.0000
## 4Hz / 33Hz 0.831 0.416 Inf 1 -0.371 1.0000
## 4Hz / 40Hz 0.804 0.408 Inf 1 -0.429 1.0000
## 4Hz / 50Hz 1.138 0.570 Inf 1 0.258 1.0000
## 4Hz / 80Hz 0.943 0.473 Inf 1 -0.116 1.0000
## 8Hz / 16Hz 1.139 0.578 Inf 1 0.257 1.0000
## 8Hz / 25Hz 1.999 1.027 Inf 1 1.348 1.0000
## 8Hz / 33Hz 1.977 1.003 Inf 1 1.343 1.0000
## 8Hz / 40Hz 1.915 0.984 Inf 1 1.264 1.0000
## 8Hz / 50Hz 2.710 1.375 Inf 1 1.964 1.0000
## 8Hz / 80Hz 2.246 1.140 Inf 1 1.594 1.0000
## 16Hz / 25Hz 1.755 0.890 Inf 1 1.108 1.0000
## 16Hz / 33Hz 1.736 0.870 Inf 1 1.101 1.0000
## 16Hz / 40Hz 1.681 0.853 Inf 1 1.024 1.0000
## 16Hz / 50Hz 2.379 1.192 Inf 1 1.730 1.0000
## 16Hz / 80Hz 1.972 0.988 Inf 1 1.355 1.0000
## 25Hz / 33Hz 0.989 0.502 Inf 1 -0.021 1.0000
## 25Hz / 40Hz 0.958 0.492 Inf 1 -0.083 1.0000
## 25Hz / 50Hz 1.356 0.688 Inf 1 0.600 1.0000
## 25Hz / 80Hz 1.124 0.570 Inf 1 0.230 1.0000
## 33Hz / 40Hz 0.969 0.492 Inf 1 -0.063 1.0000
## 33Hz / 50Hz 1.370 0.686 Inf 1 0.629 1.0000
## 33Hz / 80Hz 1.136 0.569 Inf 1 0.254 1.0000
## 40Hz / 50Hz 1.415 0.718 Inf 1 0.684 1.0000
## 40Hz / 80Hz 1.173 0.595 Inf 1 0.314 1.0000
## 50Hz / 80Hz 0.829 0.415 Inf 1 -0.375 1.0000
##
## P value adjustment: bonferroni method for 28 tests
## Tests are performed on the log scale
##
## [1] "~~~~~~~~~~~~~~ all rates ~~~~~~~~~~~~~~"
## [1] "Model summary..."
##
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.5412110414,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.91889 -1.35093 -0.49683 0.08802 1.60111
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.5754 0.3040 28.212 < 0.0000000000000002 ***
## modulation8Hz -0.7133 0.4299 -1.659 0.097070 .
## modulation16Hz -1.4230 0.4299 -3.310 0.000933 ***
## modulation25Hz -1.8539 0.4299 -4.312 0.00001616474449 ***
## modulation33Hz -2.1116 0.4300 -4.911 0.00000090540685 ***
## modulation40Hz -2.3090 0.4300 -5.370 0.00000007865917 ***
## modulation50Hz -2.5248 0.4300 -5.872 0.00000000431401 ***
## modulation80Hz -3.0032 0.4301 -6.983 0.00000000000289 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.5412) family taken to be 1)
##
## Null deviance: 281.09 on 159 degrees of freedom
## Residual deviance: 200.02 on 152 degrees of freedom
## AIC: 2473.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.5412
## Std. Err.: 0.0505
##
## 2 x log-likelihood: -2455.2150
## [1] "Model dispersion..."
## [1] 1.06478
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
##
## Response: n
## LR Chisq Df Pr(>Chisq)
## modulation 81.078 7 0.0000000000000083 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
## Estimate 2.5 % 97.5 %
## (Intercept) 5299.50000000 3082.24798764 10267.1812755
## modulation8Hz 0.49003680 0.20837906 1.1524002
## modulation16Hz 0.24099443 0.10247397 0.5667617
## modulation25Hz 0.15661855 0.06659318 0.3683466
## modulation33Hz 0.12104916 0.05146732 0.2847030
## modulation40Hz 0.09935843 0.04224335 0.2336959
## modulation50Hz 0.08007359 0.03404249 0.1883464
## modulation80Hz 0.04962732 0.02109523 0.1167501
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## modulation response SE df asymp.LCL asymp.UCL
## 4Hz 5300 1611 Inf 2921 9615
## 8Hz 2597 789 Inf 1431 4712
## 16Hz 1277 388 Inf 704 2317
## 25Hz 830 252 Inf 457 1506
## 33Hz 642 195 Inf 353 1164
## 40Hz 527 160 Inf 290 956
## 50Hz 424 129 Inf 234 770
## 80Hz 263 80 Inf 145 477
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## 4Hz / 8Hz 2.04 0.877 Inf 1 1.659 1.0000
## 4Hz / 16Hz 4.15 1.784 Inf 1 3.310 0.0261
## 4Hz / 25Hz 6.38 2.745 Inf 1 4.312 0.0005
## 4Hz / 33Hz 8.26 3.552 Inf 1 4.911 <.0001
## 4Hz / 40Hz 10.06 4.327 Inf 1 5.370 <.0001
## 4Hz / 50Hz 12.49 5.370 Inf 1 5.872 <.0001
## 4Hz / 80Hz 20.15 8.666 Inf 1 6.983 <.0001
## 8Hz / 16Hz 2.03 0.874 Inf 1 1.651 1.0000
## 8Hz / 25Hz 3.13 1.345 Inf 1 2.653 0.2233
## 8Hz / 33Hz 4.05 1.741 Inf 1 3.252 0.0321
## 8Hz / 40Hz 4.93 2.121 Inf 1 3.711 0.0058
## 8Hz / 50Hz 6.12 2.632 Inf 1 4.213 0.0007
## 8Hz / 80Hz 9.87 4.247 Inf 1 5.324 <.0001
## 16Hz / 25Hz 1.54 0.662 Inf 1 1.002 1.0000
## 16Hz / 33Hz 1.99 0.856 Inf 1 1.601 1.0000
## 16Hz / 40Hz 2.43 1.043 Inf 1 2.061 1.0000
## 16Hz / 50Hz 3.01 1.294 Inf 1 2.562 0.2912
## 16Hz / 80Hz 4.86 2.089 Inf 1 3.674 0.0067
## 25Hz / 33Hz 1.29 0.556 Inf 1 0.599 1.0000
## 25Hz / 40Hz 1.58 0.678 Inf 1 1.058 1.0000
## 25Hz / 50Hz 1.96 0.841 Inf 1 1.560 1.0000
## 25Hz / 80Hz 3.16 1.357 Inf 1 2.672 0.2112
## 33Hz / 40Hz 1.22 0.524 Inf 1 0.459 1.0000
## 33Hz / 50Hz 1.51 0.650 Inf 1 0.961 1.0000
## 33Hz / 80Hz 2.44 1.049 Inf 1 2.073 1.0000
## 40Hz / 50Hz 1.24 0.534 Inf 1 0.502 1.0000
## 40Hz / 80Hz 2.00 0.861 Inf 1 1.614 1.0000
## 50Hz / 80Hz 1.61 0.694 Inf 1 1.112 1.0000
##
## P value adjustment: bonferroni method for 28 tests
## Tests are performed on the log scale
names(all.models2) <- sort(unique(n_calls_reg$condition))
# pull out anova tables
(anova2 <- do.call(rbind, lapply(all.models2, '[[', 'deviance')) %>%
mutate(across(is.numeric, round, 2)))
## LR Chisq Df Pr(>Chisq)
## silence 0.56 7 1.00
## steady-state masker 7.34 7 0.39
## random masker 81.08 7 0.00
## modulation contrast IRR se asymp_lcl asymp_ucl z_ratio p_value
## 1 silence 4Hz / 8Hz 0.96 0.47 0.20 4.49 -0.09 1.00
## 2 silence 4Hz / 16Hz 0.91 0.45 0.19 4.24 -0.20 1.00
## 3 silence 4Hz / 25Hz 0.88 0.44 0.19 4.14 -0.25 1.00
## 4 silence 4Hz / 33Hz 0.76 0.38 0.16 3.57 -0.55 1.00
## 5 silence 4Hz / 40Hz 0.82 0.41 0.17 3.93 -0.39 1.00
## 6 silence 4Hz / 50Hz 0.96 0.48 0.21 4.50 -0.08 1.00
## 7 silence 4Hz / 80Hz 0.79 0.39 0.17 3.71 -0.47 1.00
## 8 silence 8Hz / 16Hz 0.95 0.47 0.20 4.43 -0.11 1.00
## 9 silence 8Hz / 25Hz 0.92 0.46 0.20 4.32 -0.16 1.00
## 10 silence 8Hz / 33Hz 0.80 0.39 0.17 3.72 -0.46 1.00
## 11 silence 8Hz / 40Hz 0.86 0.43 0.18 4.10 -0.30 1.00
## 12 silence 8Hz / 50Hz 1.00 0.50 0.22 4.70 0.01 1.00
## 13 silence 8Hz / 80Hz 0.83 0.41 0.18 3.87 -0.38 1.00
## 14 silence 16Hz / 25Hz 0.98 0.48 0.21 4.57 -0.05 1.00
## 15 silence 16Hz / 33Hz 0.84 0.42 0.18 3.94 -0.35 1.00
## 16 silence 16Hz / 40Hz 0.91 0.46 0.19 4.34 -0.19 1.00
## 17 silence 16Hz / 50Hz 1.06 0.52 0.23 4.97 0.12 1.00
## 18 silence 16Hz / 80Hz 0.88 0.43 0.19 4.10 -0.27 1.00
## 19 silence 25Hz / 33Hz 0.86 0.43 0.18 4.04 -0.30 1.00
## 20 silence 25Hz / 40Hz 0.93 0.47 0.20 4.44 -0.14 1.00
## 21 silence 25Hz / 50Hz 1.09 0.54 0.23 5.09 0.17 1.00
## 22 silence 25Hz / 80Hz 0.90 0.44 0.19 4.20 -0.22 1.00
## 23 silence 33Hz / 40Hz 1.08 0.54 0.23 5.16 0.15 1.00
## 24 silence 33Hz / 50Hz 1.26 0.62 0.27 5.90 0.47 1.00
## 25 silence 33Hz / 80Hz 1.04 0.51 0.22 4.87 0.08 1.00
## 26 silence 40Hz / 50Hz 1.17 0.58 0.24 5.58 0.31 1.00
## 27 silence 40Hz / 80Hz 0.96 0.48 0.20 4.60 -0.07 1.00
## 28 silence 50Hz / 80Hz 0.82 0.41 0.18 3.86 -0.39 1.00
## 29 steady-state masker 4Hz / 8Hz 0.42 0.21 0.09 2.05 -1.71 1.00
## 30 steady-state masker 4Hz / 16Hz 0.48 0.24 0.10 2.29 -1.47 1.00
## 31 steady-state masker 4Hz / 25Hz 0.84 0.43 0.17 4.10 -0.34 1.00
## 32 steady-state masker 4Hz / 33Hz 0.83 0.42 0.17 3.97 -0.37 1.00
## 33 steady-state masker 4Hz / 40Hz 0.80 0.41 0.16 3.93 -0.43 1.00
## 34 steady-state masker 4Hz / 50Hz 1.14 0.57 0.24 5.44 0.26 1.00
## 35 steady-state masker 4Hz / 80Hz 0.94 0.47 0.20 4.51 -0.12 1.00
## 36 steady-state masker 8Hz / 16Hz 1.14 0.58 0.23 5.56 0.26 1.00
## 37 steady-state masker 8Hz / 25Hz 2.00 1.03 0.40 9.95 1.35 1.00
## 38 steady-state masker 8Hz / 33Hz 1.98 1.00 0.41 9.65 1.34 1.00
## 39 steady-state masker 8Hz / 40Hz 1.92 0.98 0.38 9.54 1.26 1.00
## 40 steady-state masker 8Hz / 50Hz 2.71 1.38 0.56 13.22 1.96 1.00
## 41 steady-state masker 8Hz / 80Hz 2.25 1.14 0.46 10.96 1.59 1.00
## 42 steady-state masker 16Hz / 25Hz 1.75 0.89 0.36 8.56 1.11 1.00
## 43 steady-state masker 16Hz / 33Hz 1.74 0.87 0.36 8.30 1.10 1.00
## 44 steady-state masker 16Hz / 40Hz 1.68 0.85 0.34 8.20 1.02 1.00
## 45 steady-state masker 16Hz / 50Hz 2.38 1.19 0.50 11.37 1.73 1.00
## 46 steady-state masker 16Hz / 80Hz 1.97 0.99 0.41 9.43 1.35 1.00
## 47 steady-state masker 25Hz / 33Hz 0.99 0.50 0.20 4.83 -0.02 1.00
## 48 steady-state masker 25Hz / 40Hz 0.96 0.49 0.19 4.77 -0.08 1.00
## 49 steady-state masker 25Hz / 50Hz 1.36 0.69 0.28 6.62 0.60 1.00
## 50 steady-state masker 25Hz / 80Hz 1.12 0.57 0.23 5.48 0.23 1.00
## 51 steady-state masker 33Hz / 40Hz 0.97 0.49 0.20 4.73 -0.06 1.00
## 52 steady-state masker 33Hz / 50Hz 1.37 0.69 0.29 6.55 0.63 1.00
## 53 steady-state masker 33Hz / 80Hz 1.14 0.57 0.24 5.43 0.25 1.00
## 54 steady-state masker 40Hz / 50Hz 1.42 0.72 0.29 6.91 0.68 1.00
## 55 steady-state masker 40Hz / 80Hz 1.17 0.60 0.24 5.72 0.31 1.00
## 56 steady-state masker 50Hz / 80Hz 0.83 0.42 0.17 3.96 -0.38 1.00
## 57 random masker 4Hz / 8Hz 2.04 0.88 0.53 7.82 1.66 1.00
## 58 random masker 4Hz / 16Hz 4.15 1.78 1.08 15.89 3.31 0.03
## 59 random masker 4Hz / 25Hz 6.38 2.74 1.67 24.46 4.31 0.00
## 60 random masker 4Hz / 33Hz 8.26 3.55 2.16 31.65 4.91 0.00
## 61 random masker 4Hz / 40Hz 10.06 4.33 2.63 38.56 5.37 0.00
## 62 random masker 4Hz / 50Hz 12.49 5.37 3.26 47.85 5.87 0.00
## 63 random masker 4Hz / 80Hz 20.15 8.67 5.26 77.22 6.98 0.00
## 64 random masker 8Hz / 16Hz 2.03 0.87 0.53 7.79 1.65 1.00
## 65 random masker 8Hz / 25Hz 3.13 1.34 0.82 11.98 2.65 0.22
## 66 random masker 8Hz / 33Hz 4.05 1.74 1.06 15.51 3.25 0.03
## 67 random masker 8Hz / 40Hz 4.93 2.12 1.29 18.90 3.71 0.01
## 68 random masker 8Hz / 50Hz 6.12 2.63 1.60 23.45 4.21 0.00
## 69 random masker 8Hz / 80Hz 9.87 4.25 2.58 37.84 5.32 0.00
## 70 random masker 16Hz / 25Hz 1.54 0.66 0.40 5.89 1.00 1.00
## 71 random masker 16Hz / 33Hz 1.99 0.86 0.52 7.63 1.60 1.00
## 72 random masker 16Hz / 40Hz 2.43 1.04 0.63 9.29 2.06 1.00
## 73 random masker 16Hz / 50Hz 3.01 1.29 0.78 11.53 2.56 0.29
## 74 random masker 16Hz / 80Hz 4.86 2.09 1.27 18.61 3.67 0.01
## 75 random masker 25Hz / 33Hz 1.29 0.56 0.34 4.96 0.60 1.00
## 76 random masker 25Hz / 40Hz 1.58 0.68 0.41 6.04 1.06 1.00
## 77 random masker 25Hz / 50Hz 1.96 0.84 0.51 7.50 1.56 1.00
## 78 random masker 25Hz / 80Hz 3.16 1.36 0.82 12.10 2.67 0.21
## 79 random masker 33Hz / 40Hz 1.22 0.52 0.32 4.67 0.46 1.00
## 80 random masker 33Hz / 50Hz 1.51 0.65 0.39 5.79 0.96 1.00
## 81 random masker 33Hz / 80Hz 2.44 1.05 0.64 9.35 2.07 1.00
## 82 random masker 40Hz / 50Hz 1.24 0.53 0.32 4.76 0.50 1.00
## 83 random masker 40Hz / 80Hz 2.00 0.86 0.52 7.68 1.61 1.00
## 84 random masker 50Hz / 80Hz 1.61 0.69 0.42 6.19 1.11 1.00
# IRR
IRRs2 <- do.call(rbind, lapply(all.models2, '[[', 'IRR')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(condition = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
mutate(coef = str_remove(str_split(row.names(.),"[.]",simplify = T)[,2],"condition"), .after=1) %>%
remove_rownames()
IRRs2
## condition coef Estimate 2.5 % 97.5 %
## 1 silence (Intercept) 2981.75 1612.60 6456.97
## 2 silence modulation8Hz 1.04 0.39 2.80
## 3 silence modulation16Hz 1.10 0.41 2.96
## 4 silence modulation25Hz 1.13 0.42 3.03
## 5 silence modulation33Hz 1.31 0.49 3.52
## 6 silence modulation40Hz 1.21 0.45 3.32
## 7 silence modulation50Hz 1.04 0.39 2.79
## 8 silence modulation80Hz 1.26 0.47 3.38
## 9 steady-state masker (Intercept) 1446.90 776.51 3171.48
## 10 steady-state masker modulation8Hz 2.38 0.87 6.60
## 11 steady-state masker modulation16Hz 2.09 0.77 5.69
## 12 steady-state masker modulation25Hz 1.19 0.43 3.30
## 13 steady-state masker modulation33Hz 1.20 0.44 3.28
## 14 steady-state masker modulation40Hz 1.24 0.45 3.45
## 15 steady-state masker modulation50Hz 0.88 0.32 2.39
## 16 steady-state masker modulation80Hz 1.06 0.39 2.89
## 17 random masker (Intercept) 5299.50 3082.25 10267.18
## 18 random masker modulation8Hz 0.49 0.21 1.15
## 19 random masker modulation16Hz 0.24 0.10 0.57
## 20 random masker modulation25Hz 0.16 0.07 0.37
## 21 random masker modulation33Hz 0.12 0.05 0.28
## 22 random masker modulation40Hz 0.10 0.04 0.23
## 23 random masker modulation50Hz 0.08 0.03 0.19
## 24 random masker modulation80Hz 0.05 0.02 0.12
Use bootstrapped means from the whole data set, we now want to see whether average call onsets were “aimed” more at dodging AM peaks or more targeted at AM troughs.
To find out, we first compute two measures of call onset timing as follows:
Plot timings relative to both acoustic features. Show raw data (dots), distributions, and boxplots.
# timings by condition and modulation
more_timings <- mu_bs_time %>% group_by(modulation, condition) %>%
summarise(median_peak = median(t_from_peak),
iqr_peak = IQR(t_from_peak),
mad_peak = mad(t_from_peak),
median_trough = median(t_to_trough),
mad_trough = mad(t_to_trough),
iqr_trough = IQR(t_to_trough))
more_timings
## # A tibble: 24 × 8
## # Groups: modulation [8]
## modulation condition median_peak iqr_peak mad_peak median_trough mad_trough iqr_trough
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 4Hz silence 0.0308 0.200 0.131 0.0942 0.131 0.200
## 2 4Hz steady-state masker 0.0837 0.00310 0.00222 0.0413 0.00222 0.00310
## 3 4Hz random masker 0.0837 0.00210 0.00163 0.0413 0.00163 0.00210
## 4 8Hz silence 0.0346 0.0456 0.0249 0.0279 0.0249 0.0456
## 5 8Hz steady-state masker 0.0528 0.00130 0.000890 0.0097 0.000890 0.00130
## 6 8Hz random masker 0.0462 0.00140 0.00104 0.0163 0.00104 0.00140
## 7 16Hz silence -0.0197 0.0144 0.00875 0.0505 0.00875 0.0144
## 8 16Hz steady-state masker 0.0258 0.000900 0.000593 0.0050 0.000593 0.000900
## 9 16Hz random masker 0.0223 0.00140 0.00104 0.0085 0.00104 0.00140
## 10 25Hz silence 0.0099 0.00663 0.00482 0.0101 0.00482 0.00663
## # … with 14 more rows
pb <- p.peak_trough(mu_bs_time %>%
dplyr::filter(!condition=="silence", !condition=="half-band masker"), binwidth = 8) %>%
recolor()
print(pb)
ggsave(file.path(figures_path,paste0(exp, '_fig3a-alt1.png')), pb, width = 6, height = 4)
Run three versions of the random forest model:
A “full” model: \(modulation ~ t_{peak} + t_{trough}\)
A “peaks” model: \(modulation ~ t_{peak}\)
A “troughs” model: \(modulation ~ t_{trough}\)
mu_rf <- mu_bs_time %>% dplyr::filter(!condition=="silence") %>%
dplyr::select(t_to_trough, t_from_peak, modulation)
test_idx <- sample(1:nrow(mu_rf), size = floor(nrow(mu_rf)*0.60))
# split into 0.6:0.4 train/validation set
mu_bs_train <- mu_rf[test_idx,]
mu_bs_valid <- mu_rf[-test_idx,]
# shuffled
mu_bs_shuffled_pk <- mu_bs_train
rf_files <- c(file.path(models_path, paste0(exp,"_RF_full.RDS")),
file.path(models_path, paste0(exp,"_RF_troughs_only.RDS")),
file.path(models_path, paste0(exp,"_RF_peaks_only.RDS")))
if (all(file.exists(rf_files))) {
rf_models <- lapply(rf_files, FUN = readr::read_rds)
rf_classifier <- rf_models[[1]]
rf_classifier_tr <- rf_models[[2]]
rf_classifier_pk <- rf_models[[3]]
} else {
# train
# full model
rf_classifier <- randomForest(modulation ~ ., data=mu_bs_train, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
readr::write_rds(rf_classifier, file = rf_files[1])
# only troughs
rf_classifier_tr <- randomForest(modulation ~ t_to_trough, data=mu_bs_train, ntree=500, mtry=1, proximity=TRUE)
readr::write_rds(rf_classifier_tr, file = rf_files[2])
# only peaks
rf_classifier_pk <- randomForest(modulation ~ t_from_peak, data=mu_bs_train, ntree=500, mtry=1, proximity=TRUE)
readr::write_rds(rf_classifier_pk, file = rf_files[3])
rf_models <- list(rf_classifier, rf_classifier_tr, rf_classifier_pk)
}
rf_classifier
##
## Call:
## randomForest(formula = modulation ~ ., data = mu_bs_train, ntree = 500, mtry = 2, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 0.26%
## Confusion matrix:
## 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz 1183 0 0 0 0 0 0 0 0.000000000
## 8Hz 0 1245 0 0 0 0 0 0 0.000000000
## 16Hz 0 0 1179 0 0 0 0 0 0.000000000
## 25Hz 0 0 0 1192 15 0 0 0 0.012427506
## 33Hz 0 0 0 0 1184 2 0 0 0.001686341
## 40Hz 0 0 0 0 2 1248 3 0 0.003990423
## 50Hz 0 0 0 0 0 0 1154 3 0.002592913
## 80Hz 0 0 0 0 0 0 0 1190 0.000000000
rf_classifier_tr
##
## Call:
## randomForest(formula = modulation ~ t_to_trough, data = mu_bs_train, ntree = 500, mtry = 1, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 1
##
## OOB estimate of error rate: 36.68%
## Confusion matrix:
## 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz 1183 0 0 0 0 0 0 0 0.0000000
## 8Hz 0 848 40 16 5 80 229 27 0.3188755
## 16Hz 0 52 760 8 73 128 108 50 0.3553859
## 25Hz 0 9 99 280 86 174 342 217 0.7680199
## 33Hz 0 6 212 4 697 59 20 188 0.4123103
## 40Hz 0 1 58 1 37 652 223 281 0.4796488
## 50Hz 0 70 17 63 17 36 744 210 0.3569576
## 80Hz 0 11 86 38 79 14 47 915 0.2310924
rf_classifier_pk
##
## Call:
## randomForest(formula = modulation ~ t_from_peak, data = mu_bs_train, ntree = 500, mtry = 1, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 1
##
## OOB estimate of error rate: 15.11%
## Confusion matrix:
## 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz 1183 0 0 0 0 0 0 0 0.00000000
## 8Hz 0 1245 0 0 0 0 0 0 0.00000000
## 16Hz 0 0 1179 0 0 0 0 0 0.00000000
## 25Hz 0 0 0 730 436 5 36 0 0.39519470
## 33Hz 0 0 0 5 1129 2 2 48 0.04806071
## 40Hz 0 0 0 34 13 1010 49 147 0.19393456
## 50Hz 0 0 0 1 32 18 994 112 0.14088159
## 80Hz 0 0 0 0 12 202 297 679 0.42941176
Calculate variable importance of two predictors from the full model…
Predict modulation classes for the validation set for each of the models…
predictions <- predict(rf_classifier,mu_bs_valid[,-3])
predictionstr <- predict(rf_classifier_tr,mu_bs_valid[,-3])
predictionspk <- predict(rf_classifier_pk,mu_bs_valid[,-3])
# predicted classes vs. observed classes
table(observed=mu_bs_valid[,3],predicted=predictions)
## predicted
## observed 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 755 0 0 0 0 0 0
## 16Hz 0 0 821 0 0 0 0 0
## 25Hz 0 0 0 783 10 0 0 0
## 33Hz 0 0 0 0 811 3 0 0
## 40Hz 0 0 0 0 2 745 0 0
## 50Hz 0 0 0 0 0 2 838 3
## 80Hz 0 0 0 0 0 0 0 810
print("Full model:")
## [1] "Full model:"
# proportion misclassifications
mean(predictions != mu_bs_valid$modulation) * 100
## [1] 0.3125
# stats
(matfull <- caret::confusionMatrix(predictions, mu_bs_valid$modulation))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 755 0 0 0 0 0 0
## 16Hz 0 0 821 0 0 0 0 0
## 25Hz 0 0 0 783 0 0 0 0
## 33Hz 0 0 0 10 811 2 0 0
## 40Hz 0 0 0 0 3 745 2 0
## 50Hz 0 0 0 0 0 0 838 0
## 80Hz 0 0 0 0 0 0 3 810
##
## Overall Statistics
##
## Accuracy : 0.9969
## 95% CI : (0.9952, 0.9981)
## No Information Rate : 0.1317
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.9964
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity 1.0000 1.000 1.0000 0.9874 0.9963 0.9973 0.9941
## Specificity 1.0000 1.000 1.0000 1.0000 0.9979 0.9991 1.0000
## Pos Pred Value 1.0000 1.000 1.0000 1.0000 0.9854 0.9933 1.0000
## Neg Pred Value 1.0000 1.000 1.0000 0.9982 0.9995 0.9996 0.9991
## Prevalence 0.1277 0.118 0.1283 0.1239 0.1272 0.1167 0.1317
## Detection Rate 0.1277 0.118 0.1283 0.1223 0.1267 0.1164 0.1309
## Detection Prevalence 0.1277 0.118 0.1283 0.1223 0.1286 0.1172 0.1309
## Balanced Accuracy 1.0000 1.000 1.0000 0.9937 0.9971 0.9982 0.9970
## Class: 80Hz
## Sensitivity 1.0000
## Specificity 0.9995
## Pos Pred Value 0.9963
## Neg Pred Value 1.0000
## Prevalence 0.1266
## Detection Rate 0.1266
## Detection Prevalence 0.1270
## Balanced Accuracy 0.9997
print("Troughs model:")
## [1] "Troughs model:"
# stats
(mattr <- caret::confusionMatrix(predictionstr, mu_bs_valid$modulation))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 520 35 6 4 0 37 12
## 16Hz 0 23 539 64 156 39 14 41
## 25Hz 0 13 3 193 4 1 52 39
## 33Hz 0 3 56 50 449 9 15 65
## 40Hz 0 49 91 104 49 390 26 5
## 50Hz 0 123 67 215 17 130 568 41
## 80Hz 0 24 30 161 135 178 131 607
##
## Overall Statistics
##
## Accuracy : 0.638
## 95% CI : (0.6261, 0.6498)
## No Information Rate : 0.1317
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.5858
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity 1.0000 0.68874 0.65652 0.24338 0.55160 0.52209 0.67378
## Specificity 1.0000 0.98335 0.93959 0.98002 0.96455 0.94269 0.89329
## Pos Pred Value 1.0000 0.84691 0.61530 0.63279 0.69397 0.54622 0.48923
## Neg Pred Value 1.0000 0.95938 0.94895 0.90156 0.93655 0.93721 0.94751
## Prevalence 0.1277 0.11797 0.12828 0.12391 0.12719 0.11672 0.13172
## Detection Rate 0.1277 0.08125 0.08422 0.03016 0.07016 0.06094 0.08875
## Detection Prevalence 0.1277 0.09594 0.13687 0.04766 0.10109 0.11156 0.18141
## Balanced Accuracy 1.0000 0.83604 0.79806 0.61170 0.75808 0.73239 0.78354
## Class: 80Hz
## Sensitivity 0.74938
## Specificity 0.88211
## Pos Pred Value 0.47946
## Neg Pred Value 0.96046
## Prevalence 0.12656
## Detection Rate 0.09484
## Detection Prevalence 0.19781
## Balanced Accuracy 0.81575
print("Peaks model:")
## [1] "Peaks model:"
# stats
(matpk <- caret::confusionMatrix(predictionspk, mu_bs_valid$modulation))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 755 0 0 0 0 0 0
## 16Hz 0 0 821 0 0 0 0 0
## 25Hz 0 0 0 471 4 26 0 0
## 33Hz 0 0 0 290 771 8 19 0
## 40Hz 0 0 0 8 0 603 17 150
## 50Hz 0 0 0 23 4 33 747 190
## 80Hz 0 0 0 1 35 77 60 470
##
## Overall Statistics
##
## Accuracy : 0.8523
## 95% CI : (0.8434, 0.861)
## No Information Rate : 0.1317
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.8312
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity 1.0000 1.000 1.0000 0.59395 0.9472 0.80723 0.8861
## Specificity 1.0000 1.000 1.0000 0.99465 0.9433 0.96904 0.9550
## Pos Pred Value 1.0000 1.000 1.0000 0.94012 0.7086 0.77506 0.7492
## Neg Pred Value 1.0000 1.000 1.0000 0.94541 0.9919 0.97439 0.9822
## Prevalence 0.1277 0.118 0.1283 0.12391 0.1272 0.11672 0.1317
## Detection Rate 0.1277 0.118 0.1283 0.07359 0.1205 0.09422 0.1167
## Detection Prevalence 0.1277 0.118 0.1283 0.07828 0.1700 0.12156 0.1558
## Balanced Accuracy 1.0000 1.000 1.0000 0.79430 0.9452 0.88814 0.9206
## Class: 80Hz
## Sensitivity 0.58025
## Specificity 0.96905
## Pos Pred Value 0.73095
## Neg Pred Value 0.94094
## Prevalence 0.12656
## Detection Rate 0.07344
## Detection Prevalence 0.10047
## Balanced Accuracy 0.77465
validation <- list(matfull, mattr, matpk)
# training matrices
# confmats <- lapply(rf_models,function(x) cormat_to_df(x$confusion[,-3])) %>%
# setNames(c("full","troughs","peaks")) %>%
# map_df(., ~as.data.frame(.x), .id="model") %>%
# mutate(model = factor(model, levels = c("troughs","peaks","full")))
# validation matrices
confmats <- lapply(validation, function(x) cormat_to_df(x$table %>% as.data.frame()) ) %>%
setNames(c("full", "troughs", "peaks")) %>%
map_df(., ~as.data.frame(.x), .id="model") %>%
mutate(model = factor(model, levels = c("troughs", "peaks", "full")))
confmat <- p.mat(confmats, ulim = max((mu_bs_valid%>%group_by(modulation)%>%summarise(n=n()))[,2])) +
facet_grid(cols=vars(model))
print(confmat)
ggsave(file.path(figures_path,paste0(exp, '_fig3b.png')),
confmat, width = h, height = w2)
Now, let’s see what happens if we shuffle the modulation class labelled for one or the other measure in our full model. The drop in performance, respectively, should give us some insight into which variable is doing more of the legwork.
# shuffled
mu_bs_shuffled_pk <- mu_bs_train
mu_bs_shuffled_pk$t_from_peak <- mu_bs_shuffled_pk$t_from_peak[sample(1:nrow(mu_bs_shuffled_pk))]
mu_bs_shuffled_tr <- mu_bs_train
mu_bs_shuffled_tr$t_to_trough <- mu_bs_shuffled_tr$t_to_trough[sample(1:nrow(mu_bs_shuffled_tr))]
# train
# full model
shPeak_rf_classifier <- randomForest(modulation ~ ., data=mu_bs_shuffled_pk, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
shPeak_rf_classifier
##
## Call:
## randomForest(formula = modulation ~ ., data = mu_bs_shuffled_pk, ntree = 500, mtry = 2, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 42.24%
## Confusion matrix:
## 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz 1183 0 0 0 0 0 0 0 0.0000000
## 8Hz 0 865 58 21 13 57 197 34 0.3052209
## 16Hz 0 60 569 91 171 82 108 98 0.5173876
## 25Hz 0 19 98 299 76 112 328 275 0.7522784
## 33Hz 0 9 162 81 601 33 32 268 0.4932546
## 40Hz 0 5 93 116 37 489 210 303 0.6097366
## 50Hz 0 105 17 69 30 42 655 239 0.4338807
## 80Hz 0 8 49 48 131 13 57 884 0.2571429
shTrough_rf_classifier <- randomForest(modulation ~ ., data=mu_bs_shuffled_tr, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
shTrough_rf_classifier
##
## Call:
## randomForest(formula = modulation ~ ., data = mu_bs_shuffled_tr, ntree = 500, mtry = 2, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 17.32%
## Confusion matrix:
## 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz 1183 0 0 0 0 0 0 0 0.0000000
## 8Hz 0 1245 0 0 0 0 0 0 0.0000000
## 16Hz 0 0 1179 0 0 0 0 0 0.0000000
## 25Hz 0 0 0 757 411 30 9 0 0.3728252
## 33Hz 0 0 0 18 1095 6 29 38 0.0767285
## 40Hz 0 0 0 26 14 962 52 199 0.2322426
## 50Hz 0 0 0 1 34 35 882 205 0.2376837
## 80Hz 0 0 0 0 33 224 299 634 0.4672269
Let’s get predictions as before…
predictionstr_sh <- predict(shTrough_rf_classifier,mu_bs_valid[,-3])
predictionspk_sh <- predict(shPeak_rf_classifier,mu_bs_valid[,-3])
print("Shuffled troughs:")
## [1] "Shuffled troughs:"
# stats
shtr <- caret::confusionMatrix(predictionstr_sh , mu_bs_valid$modulation)
shtr
## Confusion Matrix and Statistics
##
## Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 755 0 0 0 0 0 0
## 16Hz 0 0 821 0 0 0 0 0
## 25Hz 0 0 0 502 24 24 1 0
## 33Hz 0 0 0 264 744 7 19 28
## 40Hz 0 0 0 8 2 616 27 111
## 50Hz 0 0 0 18 22 23 637 178
## 80Hz 0 0 0 1 22 77 159 493
##
## Overall Statistics
##
## Accuracy : 0.8414
## 95% CI : (0.8322, 0.8503)
## No Information Rate : 0.1317
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.8187
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity 1.0000 1.000 1.0000 0.63304 0.9140 0.82463 0.75563
## Specificity 1.0000 1.000 1.0000 0.99126 0.9431 0.97382 0.95663
## Pos Pred Value 1.0000 1.000 1.0000 0.91107 0.7006 0.80628 0.72551
## Neg Pred Value 1.0000 1.000 1.0000 0.95025 0.9869 0.97676 0.96269
## Prevalence 0.1277 0.118 0.1283 0.12391 0.1272 0.11672 0.13172
## Detection Rate 0.1277 0.118 0.1283 0.07844 0.1163 0.09625 0.09953
## Detection Prevalence 0.1277 0.118 0.1283 0.08609 0.1659 0.11937 0.13719
## Balanced Accuracy 1.0000 1.000 1.0000 0.81215 0.9285 0.89923 0.85613
## Class: 80Hz
## Sensitivity 0.60864
## Specificity 0.95367
## Pos Pred Value 0.65559
## Neg Pred Value 0.94387
## Prevalence 0.12656
## Detection Rate 0.07703
## Detection Prevalence 0.11750
## Balanced Accuracy 0.78115
print("Shuffled peaks:")
## [1] "Shuffled peaks:"
# stats
shpk <- caret::confusionMatrix(predictionspk_sh , mu_bs_valid$modulation)
shpk
## Confusion Matrix and Statistics
##
## Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
## 4Hz 817 0 0 0 0 0 0 0
## 8Hz 0 459 37 0 4 0 88 0
## 16Hz 0 57 412 61 107 59 6 19
## 25Hz 0 25 94 168 60 18 50 40
## 33Hz 0 5 73 59 394 14 29 56
## 40Hz 0 48 66 81 26 357 28 13
## 50Hz 0 134 74 236 33 95 494 57
## 80Hz 0 27 65 188 190 204 148 625
##
## Overall Statistics
##
## Accuracy : 0.5822
## 95% CI : (0.57, 0.5943)
## No Information Rate : 0.1317
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.522
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity 1.0000 0.60795 0.50183 0.21185 0.48403 0.47791 0.58600
## Specificity 1.0000 0.97715 0.94461 0.94881 0.95775 0.95365 0.88681
## Pos Pred Value 1.0000 0.78061 0.57143 0.36923 0.62540 0.57674 0.43989
## Neg Pred Value 1.0000 0.94907 0.92798 0.89487 0.92721 0.93254 0.93386
## Prevalence 0.1277 0.11797 0.12828 0.12391 0.12719 0.11672 0.13172
## Detection Rate 0.1277 0.07172 0.06438 0.02625 0.06156 0.05578 0.07719
## Detection Prevalence 0.1277 0.09187 0.11266 0.07109 0.09844 0.09672 0.17547
## Balanced Accuracy 1.0000 0.79255 0.72322 0.58033 0.72089 0.71578 0.73641
## Class: 80Hz
## Sensitivity 0.77160
## Specificity 0.85295
## Pos Pred Value 0.43193
## Neg Pred Value 0.96265
## Prevalence 0.12656
## Detection Rate 0.09766
## Detection Prevalence 0.22609
## Balanced Accuracy 0.81228
validation_sh <- list(shtr, shpk)
How did variable importance metrics compare between the full (unshuffled) version and the two shuffled versions?
importance_sh <- importance(shTrough_rf_classifier) %>% importance_to_df()
importance_shpk <- importance(shPeak_rf_classifier) %>% importance_to_df()
compare <- rbind(importance_shpk %>% mutate(shuffle = "shuffled peaks"),
importance_sh %>% mutate(shuffle = "shuffled troughs"),
importance %>% mutate(shuffle = "original"))
p.imp.sh <- p.importance_shuffled(compare)
p.imp.sh
ggsave(file.path(figures_path,paste0(exp, '_fig3c.png')),p.imp.sh, width = 4, height = 4)
Finally, here’s how our shuffled modelled fared at classifying modulation rate classes…
# validation matrices
confmats_sh <- lapply(validation_sh, function(x) cormat_to_df(x$table %>% as.data.frame()) ) %>%
setNames(c("troughs", "peaks")) %>%
map_df(., ~as.data.frame(.x), .id="model") %>%
mutate(model = factor(model, levels = c("troughs", "peaks", "full")))
confmat_sh <- p.mat(confmats_sh, ulim = max(confmats_sh$n)) +
facet_grid(cols=vars(model))
print(confmat_sh)
ggsave(file.path(figures_path,paste0(exp, '_figs6.png')),
confmat_sh, width = h, height = w3)
overlap_calls <- data[which(data$onset_interval<data$duration),]%>%
mutate(overlap = duration-onset_interval)
cat("Of a total of", nrow(data), "calls,", nrow(overlap_calls),
"(", round(nrow(overlap_calls)/nrow(data)*100,2), "%) overlapped in time.")
## Of a total of 1088994 calls, 31942 ( 2.93 %) overlapped in time.
n_overlap <- overlap_calls %>% group_by(modulation, condition) %>% summarise(n = n())
(n_overlap_conttable <- addmargins(table(overlap_calls$modulation,overlap_calls$condition), c(1,2)))
##
## silence steady-state masker random masker Sum
## 4Hz 3514 406 1331 5251
## 8Hz 3461 1322 603 5386
## 16Hz 2452 847 330 3629
## 25Hz 3278 383 221 3882
## 33Hz 3506 345 164 4015
## 40Hz 2988 338 130 3456
## 50Hz 3017 219 99 3335
## 80Hz 2687 233 68 2988
## Sum 24903 4093 2946 31942
Let’s run the overlap data through a negative binomial…
Deviance tables…
(ovlp.anova <- do.call(rbind, lapply(ovlp.models, '[[', 'deviance')) %>%
mutate(across(is.numeric, round, 2)))
## LR Chisq Df Pr(>Chisq)
## condition 5.23 2 0.07
## condition1 7.27 2 0.03
## condition2 12.19 2 0.00
## condition3 15.20 2 0.00
## condition4 25.20 2 0.00
## condition5 31.39 2 0.00
## condition6 34.76 2 0.00
## condition7 35.55 2 0.00
And contrasts…
(contrasts <- do.call(rbind, lapply(ovlp.models, '[[', 'contrasts')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(modulation = factor(rep(unique(n_overlap$modulation),each=length(unique(n_overlap$condition)))),.before=1) %>%
remove_rownames() %>%
dplyr::select(-c("df", "null")) %>%
dplyr::rename(IRR = ratio))
## modulation masking_condition IRR se asymp_lcl asymp_ucl z_ratio p_value
## 1 4Hz silence / (steady-state masker) 4.58 3.04 0.94 22.45 2.29 0.07
## 2 4Hz silence / random masker 2.33 1.33 0.60 9.11 1.49 0.41
## 3 4Hz (steady-state masker) / random masker 0.51 0.34 0.10 2.58 -1.00 0.96
## 4 8Hz silence / (steady-state masker) 1.92 1.21 0.43 8.65 1.04 0.90
## 5 8Hz silence / random masker 5.36 3.16 1.31 21.96 2.85 0.01
## 6 8Hz (steady-state masker) / random masker 2.79 1.78 0.60 12.88 1.61 0.32
## 7 16Hz silence / (steady-state masker) 2.21 1.16 0.63 7.77 1.52 0.39
## 8 16Hz silence / random masker 6.56 3.32 1.95 22.02 3.72 0.00
## 9 16Hz (steady-state masker) / random masker 2.96 1.60 0.81 10.83 2.00 0.14
## 10 25Hz silence / (steady-state masker) 5.14 3.10 1.21 21.80 2.71 0.02
## 11 25Hz silence / random masker 9.89 5.80 2.43 40.26 3.91 0.00
## 12 25Hz (steady-state masker) / random masker 1.93 1.27 0.40 9.37 0.99 0.96
## 13 33Hz silence / (steady-state masker) 5.72 3.07 1.58 20.69 3.24 0.00
## 14 33Hz silence / random masker 14.70 7.47 4.35 49.66 5.28 0.00
## 15 33Hz (steady-state masker) / random masker 2.57 1.50 0.64 10.41 1.62 0.32
## 16 40Hz silence / (steady-state masker) 6.48 3.29 1.92 21.87 3.68 0.00
## 17 40Hz silence / random masker 18.39 9.20 5.55 60.93 5.82 0.00
## 18 40Hz (steady-state masker) / random masker 2.84 1.53 0.78 10.35 1.93 0.16
## 19 50Hz silence / (steady-state masker) 11.27 6.43 2.88 44.16 4.25 0.00
## 20 50Hz silence / random masker 27.70 15.51 7.26 105.80 5.93 0.00
## 21 50Hz (steady-state masker) / random masker 2.46 1.45 0.60 10.12 1.52 0.38
## 22 80Hz silence / (steady-state masker) 7.69 3.88 2.30 25.73 4.04 0.00
## 23 80Hz silence / random masker 26.34 13.56 7.68 90.35 6.35 0.00
## 24 80Hz (steady-state masker) / random masker 3.43 1.94 0.89 13.24 2.18 0.09
And finally, predict values…
ovlp_preds <- do.call(rbind, lapply(ovlp.models, '[[', 'predictions')) %>%
mutate(across(is.numeric, round, 2)) %>%
mutate(modulation = factor(rep(unique(n_overlap$modulation),each=length(unique(n_overlap$condition))))) %>%
remove_rownames()
ovlp_pred <- p.n_pred(ovlp_preds) %>% recolor()
ovlp_pred
ggsave(file.path(figures_path,paste0(exp, '_figs7b.png') ), ovlp_pred, width = 5, height = 4)
Finally, we can calculate circular stats for overlapping calls, and visualize the angular vectors representing those distributions for each condition:
overlap_circsummary <- overlap_calls %>%
group_by(modulation, condition) %>%
summarize(n = n(),
theta_bar = mean.circular(start_phase_circ, na.rm = T),
r_bar = rho.circular(start_phase_circ, na.rm = T),
vm = var.circular(start_phase_circ),
rt = rayleigh.test(start_phase_circ, mu=circular(0))$statistic,
p = rayleigh.test(start_phase_circ, mu=circular(0))$p.value,
bs_mu_ci_l = mle.vonmises.bootstrap.ci(start_phase_circ, alpha=.05)$mu.ci[1],
bs_mu_ci_h = mle.vonmises.bootstrap.ci(start_phase_circ, alpha=.05)$mu.ci[2],) %>%
mutate(p_adjust = p.adjust(p, method ="bonferroni"), .after = "p") %>%
mutate(across(where(is.numeric), round, 3))
overlap_circsummary %>%
kbl(escape = FALSE,
booktabs = TRUE,
format = "html",
caption = "Circular statistics for overlapping calls") %>%
kable_styling() %>%
cat(., file = file.path(tables_path,paste0(exp, "_circ_stats_overlap.html")))
p.ang <- p.angular_vectors(overlap_circsummary %>% confint_wrap())
p.ang
ggsave(file.path(figures_path,paste0(exp, '_figs7a.png') ), p.ang, width = w, height = 4)
po <- p.circ_density(dat=overlap_calls, circ_dat=overlap_circsummary %>% confint_wrap(), bins = 30) %>% recolor(which = "f")
print(po)
ggsave(file.path(figures_path,paste0(exp, '_figs7a-alt2.png') ), po, width = w, height = 5)
To determine how many cycles of continuous stimulation are needed for the bats’s timing adaptation to stabilize, estimate mean and concentration parameters for successively larger cycle windows:
max_n <- data_init %>% group_by(modulation, condition) %>% summarise(n = max(start_cycle)) %>%
ungroup() %>% summarise(min = min(n))
by_cycle_summary <- data.frame()
for (c in seq(0, max_n$min, 10)) {
#for (c in 2:length(seq(0, max_n, 10))) {
#s <- seq(0, 1000, 10)
dat_summary <- data_init %>%
dplyr::filter(start_cycle <= c) %>%
#dplyr::filter(start_cycle > s[c-1] & start_cycle <= s[c]) %>% # non-cumulative
group_by(modulation, condition) %>% # pool over groups
summarize(n = c,
theta_bar = mean.circular(start_phase_circ, na.rm = T) %>% as.numeric(),
r_bar = rho.circular(start_phase_circ, na.rm = T))
by_cycle_summary <- rbind(by_cycle_summary, dat_summary)
}
Some references for circular stats can be found from the NCSS website, as well as the documentation for the Matlab library for circular stats, and articles written on Bayesian modelling of circular data.
There is also the textbook, Statistical analysis of circular data, as well as: - Batschelet, E. (1981). Circular Statistics in Biology. Academic Press: London. - Jammalamadaka, S.R. & Sengupta, A. (2001). Topics in Circular Statistics. World Scientific, River Edge, N.J.